PGF90 (Version 12.8) 08/30/2020 00:27:52 page 1 Switches: -noasm -nodclchk -nodebug -nodlines -noline -list -idir ../include -inform warn -opt 1 -nosave -object -noonetrip -depchk on -nostandard -nosymbol -noupcase Filename: slp.f90 ( 1) subroutine slp ( 2) ( 3) USE mod_vrbls ( 4) USE mod_extra ( 5) USE mod_masks ( 6) implicit none ( 7) include 'param_o.h' ( 8) !!JNT include 'extra.comm' ../include/param_o.h ( 1)* integer,parameter :: im0=401 ( 2)* integer, parameter :: nm=6 ( 3)* integer, parameter :: lm = 50 ( 4)* integer, parameter :: nsub = 10 ( 5)* ( 6)* integer, parameter :: im=im0, jm=im ( 7)* integer, parameter :: im1=im-1, jm1=jm-1 ( 8)* ( 9)* integer, parameter :: lm1 = lm-1, lp1 =lm +1 ( 10)* ( 11)* integer, parameter :: ixm = nsub, jym = ixm, nxy = ixm*jym ( 12)* integer, parameter :: ildom = (im - 1)/ixm, jldom = (jm - 1)/jym ( 13)* integer, parameter :: ilm = (im - 1)/ixm +1, jlm = (jm - 1)/jym +1 ( 14)* ( 15)* logical,parameter::flat=.false. ( 16)* logical,parameter::hstst=.false. ( 17)* ( 18)*! integer,parameter::igm=360, jgm = 181 ( 19)* integer,parameter::igm=(im0-1)*4, jgm = (igm/2)+1 ( 20)* real,parameter::alfa=0.0, beta=0.0, gamm=0.0 ( 21)*! real,parameter::alfa=0.*3.1415926/180.,beta=66.*3.1415926/180. & ( 22)*! ,gamm=175.*3.1415926/180. ( 23)* ( 24)*!GSM integer,parameter::lsm=20 ( 25)* integer,parameter::lsm=31 ( 26)* ( 27)*!RESTART ( 28)* ( 29)* character(len=10):: restartdate='2020080412' !restart date to create appropriate folder name ( 9) include 'mapot.comm' ../include/mapot.comm ( 1)* real,dimension(lsm)::spl,alsl ( 2)* ( 3)* common/mapot/ spl,alsl ( 10) include 'const.h' ../include/const.h ( 1)* logical :: run, first, restrt, subpost ( 2)* integer :: nfcst,nbc,list,ntsd,nddamp,nprec, & ( 3)* nboco,nshde,ncp,ntddmp ( 4)* logical,parameter::lcornerm=.FALSE. ( 5)* ( 6)* real,parameter :: dt=40 ! length of each time step PGF90 (Version 12.8) 08/30/2020 00:27:52 page 2 ( 7)* logical,parameter:: sigma=.false. ( 8)*!GSM integer,parameter::outnum=40 ! number of outputs ( 9)*!GSM integer,parameter::nday=10 ! number of days ( 10)* integer,parameter:: ntstm=3600*24*10/dt ! total time steps ( 11)*! integer,parameter:: idtad=1 ( 12)*! integer,parameter:: ncnvc=45 ( 13)*!test integer,parameter::nphs=45 ( 14)* ( 15)* integer,parameter:: idtad=2 ( 16)* integer,parameter:: ncnvc=6 ( 17)* integer,parameter::nphs=6 ( 18)* ( 19)* ( 20)* integer,parameter::nradsh=1 ( 21)*!GSM integer,parameter::nradlh=2 ( 22)* integer,parameter::nradlh=1 ( 23)* real,parameter::tsph=3600./dt ( 24)* integer,parameter::nrads=tsph*nradsh ( 25)* integer,parameter::nradl=tsph*nradlh ( 26)* real,parameter::weig=0.25 ( 27)* ( 28)*! ( 29)*! parameters for diagnostics ( 30)*! ( 31)* ( 32)* integer,parameter::idgns=10 ( 33)* integer,parameter::dgnstr=24*3600*200/dt ( 34)* integer,parameter::dgnnum=(ntstm-dgnstr)/idgns+1 ( 35)* integer,parameter::ickmm=1 ( 36)* ( 37)*! ( 38)*! parameter for initial values ( 39)*! ( 40)* real(KIND=4),parameter::pt = 2500.0 ( 41)*!!!test Dragan 1dec2015 ( 42)*!! real,parameter::pt=1000 ( 43)* ( 44)*! ( 45)*! date ( 46)*! ( 47)* integer::idat(3) ( 48)*!GSM data idat/07,21,2020/ ( 49)*!GSM integer::ihrst=0 ( 50)* ( 51)*! variable sst ( 52)* logical::lvsst ( 53)* ( 54)*! data assimilation constants ( 55)* ( 56)* integer,parameter::ndassim=3*3600/dt ( 57)*! integer,parameter::ndassimm=8*ndassim ( 58)* integer,parameter::ndassimm=0 ( 59)* integer,parameter::outstr=ndassimm ! step at which starts output ( 60)* integer,parameter::outend= ntstm+ndassimm ! step at which ends output ( 61)*!GSM integer,parameter::iout=(outend-outstr)/outnum ! steps between two output ( 62)* ( 63)*! data choice ( 64)* PGF90 (Version 12.8) 08/30/2020 00:27:52 page 3 ( 65)* integer,parameter::sstc=1 ! 1 is NCEP sst data ( 66)* ! 2 is TMI sst data ( 67)* ( 68)* ( 11) include 'dynam.comm' ../include/dynam.comm ( 1)* real::fadv,fadv2, fadt, rd, f4d, ef4t, fkin ( 2)* real, dimension (lm) :: deta, rdeta, aeta, daeta,f4q2 ( 3)* real, dimension (lp1) :: eta, dfl ( 4)* ( 5)* common/dynam/ fadv,fadv2,fadt,rd,f4d,f4q2,ef4t,fkin,deta,rdeta, & ( 6)* eta,dfl,aeta,daeta ( 7)* ( 8)* ( 9)*!!!!!!!!!!!copy from dynam_comm.h ( 10)*! real :: fadv, fadv2,fadt, rd, f4d, ef4t, fkin, fcp ( 11)*! real, dimension (lm) :: deta, rdeta, aeta, daeta,f4q2 ( 12)*! real, dimension (lp1) :: eta, dfl ( 13)*! real, dimension (0:im+1,0:jm+1,nm) :: wpdar, f11, f12, f21, f22, & ( 14)*! p11, p12, p21, p22, fdiv, fddmp,fvdiff & ( 15)*! ,hbmsk,hsinp,hcosp ( 16)*! common /dynam/ fadv,fadv2, fadt, rd, f4d,f4q2, ef4t, fkin, & ( 17)*! deta, rdeta, aeta, eta, dfl, daeta, & ( 18)*! wpdar, f11, f12, f21, f22, & ( 19)*! p11, p12, p21, p22, fcp, fdiv, fddmp,hsinp,hcosp,& ( 20)*! fvdiff,hbmsk ( 12) include 'params' ( 13) ../include/params ( 1)* real,parameter::p1000=1000.e2,CAPA=0.28589641E0 ( 14) integer,parameter::nrlx1=500 ( 15) integer::i,j,n,l,lhmnt,lmap1,ll,nrlx,kmn,kmntm,lmst,li,kmm, & ( 16) km,n1 ( 17) integer::isCorner2,lma ( 18) real,dimension(0:im+1,0:jm+1,nm)::ttv,pdsl1,pbi,slpx ( 19) integer,dimension((im+2)*(jm+2)*nm)::imnt,jmnt,nmnt ( 20) real::tgss,pbin,phbi,ptin,trtv,phti,dposp,prbin,prtin,alpp1,slop, & ( 21) slpp,ttt ( 22) real,parameter::OVERRC=1.50,AD05=OVERRC*0.05,CFT0=OVERRC-1., & ( 23) r=287.04,rog=r/9.8 ( 24) logical:: stdrd ( 25) ( 26) if(sigma)then ( 27) stdrd=.true. ( 28) else ( 29) stdrd=.false. ( 30) endif ( 31) ! ( 32) !----------------------------------------------------------------------- ( 33) ! ( 34) !*** FIND THE MOST ELEVATED GLOBAL LAYER WHERE THERE IS LAND ( 35) ! ( 36) !----------------------------------------------------------------------- ( 37) DO L=1,LM ( 38) do n=1,nm ( 39) DO J=1,jm ( 40) DO I=1,IM ( 41) IF(htm(I,J,n,l).LT.0.5)GO TO 666 PGF90 (Version 12.8) 08/30/2020 00:27:52 page 4 ( 42) ENDDO ( 43) ENDDO ( 44) ENDDO ( 45) ! ( 46) !*** IF WE GET TO HERE, WE HAVE ALL ATMOSPHERE ( 47) ! ( 48) GOTO 667 ( 49) 666 CONTINUE ( 50) ! ( 51) !*** IF WE GET HERE, WE HAVE FOUND A NON_ATM POINT ( 52) ! ( 53) LHMNT=L ( 54) GOTO 668 ( 55) 667 CONTINUE ( 56) ENDDO ( 57) ! ( 58) LHMNT=LM+1 ( 59) !cccccc GO TO 669 ( 60) 668 CONTINUE ( 61) ( 62) !----------------------------------------------------------------------- ( 63) !*** ( 64) !*** INITIALIZE ARRAYS. LOAD SLP ARRAY WITH SURFACE PRESSURE. ( 65) !*** ( 66) do n=1,nm ( 67) DO J=1,jm ( 68) DO I=1,IM ( 69) PSLP(I,J,n)=0. ( 70) TTV(I,J,n)=0. ( 71) ENDDO ( 72) ENDDO ( 73) ENDDO ( 74) ! ( 75) do 110 n=1,nm ( 76) DO 110 J=1,jm ( 77) DO 110 I=1,IM ( 78) PDSL1(I,J,n)=RES(I,J,n)*PD(I,J,n) ( 79) PSLP(I,J,n)=PD(I,J,n)+PT ( 80) PBI (I,J,n)=PSLP(I,J,n) ( 81) 110 CONTINUE ( 82) ( 83) if(.not.stdrd)then ( 84) ( 85) LL=LM ( 86) ! ( 87) do n=1,nm ( 88) DO J=1,jm ( 89) DO I=1,im ( 90) IF(HTM(I,J,n,LL).LE.0.5) THEN ( 91) TGSS=FIS(I,J,n)/(R*ALOG((PDSL1(I,J,n)+PT)/(PD(I,J,n)+PT))) ( 92) LMAP1=LMH(I,J,n)+1 ( 93) ! ( 94) if(n.eq.3.and.i.eq.30.and.j.eq.101)then ( 95) ! print *,i,j,tgss,fis(i,j,n),r,pdsl1(i,j,n),pd(i,j,n) ( 96) ! stop ( 97) endif ( 98) DO 260 L=LMAP1,LM ( 99) T(I,J,n,L)=TGSS PGF90 (Version 12.8) 08/30/2020 00:27:52 page 5 ( 100) 260 CONTINUE ( 101) ! ( 102) ENDIF ( 103) ENDDO ( 104) ENDDO ( 105) ENDDO ( 106) ( 107) !---------------------------------------------------------------- ( 108) ! ( 109) !*** CREATE A TEMPORARY TV ARRAY, AND FOLLOW BY SEQUENTIAL ( 110) !*** OVERRELAXATION, DOING NRLX PASSES. ( 111) ! ( 112) !---------------------------------------------------------------- ( 113) NRLX=NRLX1 ( 114) ! ( 115) !---------------------------------------------------------------- ( 116) !---------------------------------------------------------------- ( 117) DO 300 L=LHMNT,LM ( 118) !---------------------------------------------------------------- ( 119) !---------------------------------------------------------------- ( 120) ! ( 121) KMN=0 ( 122) KMNTM=0 ( 123) ! ( 124) do 240 n=1,nm ( 125) DO 240 J=1,jm ( 126) DO 240 I=1,im ( 127) IF(HTM(I,J,n,L).GT.0.5)GO TO 240 ( 128) KMN=KMN+1 ( 129) IMNT(KMN)=I ( 130) JMNT(KMN)=J ( 131) nmnt(kmn)=n ( 132) 240 CONTINUE ( 133) ! ( 134) KMNTM=KMN ( 135) ! ( 136) do 270 n=1,nm ( 137) DO 270 J=1,jm ( 138) DO 270 I=1,IM ( 139) TTV(I,J,n)=T(I,J,n,L) ( 140) 270 CONTINUE ( 141) ! ( 142) !---------------------------------------------------------------- ( 143) !*** FOR GRID BOXES NEXT TO MOUNTAINS REPLACE TTV BY AN "EQUIVALENT" ( 144) !*** TV, ONE WHICH CORRESPONDS TO THE CHANGE IN P BETWEEN REFERENCE ( 145) !*** INTERFACE GEOPOTENTIALS, INSTEAD OF BETWEEN LAYER INTERFACES ( 146) !---------------------------------------------------------------- ( 147) ! ( 148) do n=1,nm ( 149) DO J=1,jm ( 150) DO I=1,im ( 151) IF(HTM(I,J,n,L).GT.0.5.AND. & ( 152) HTM(I,J-1,n,L)*HTM(I+1,J-1,n,L) & ( 153) *HTM(I,J+1,n,L)*HTM(I+1,J+1,n,L) & ( 154) *HTM(I-1 ,J ,n,L)*HTM(I+1 ,J ,n,L) & ( 155) *HTM(I-1 ,J-1,n,L)*HTM(I-1 ,J+1,n,L).LT.0.5)THEN ( 156) LMST=LMH(I,J,n) ( 157) !*** PGF90 (Version 12.8) 08/30/2020 00:27:52 page 6 ( 158) !*** FIND P AT THE REFERENCE INTERFACE GEOPOTENTIAL AT THE BOTTOM ( 159) !*** ( 160) PBIN=PT+PD(I,J,n) ( 161) PHBI=DFL(LMST+1) ( 162) ! ( 163) DO LI=LMST,1,-1 ( 164) PTIN=PBIN-DETA(LI)*PD(I,J,n)*RES(I,J,n) ( 165) TRTV=2.*R*T(I,J,n,LI)*(1.+0.608*Q(I,J,n,LI)) ( 166) PHTI=PHBI+TRTV*(PBIN-PTIN)/(PBIN+PTIN) ( 167) if(i.eq.17.and.j.eq.89.and.n.eq.5)then ( 168) ! print *,l,pbin,ptin,deta(li),res(i,j,n),pd(i,j,n) ( 169) endif ( 170) IF(PHTI.GE.DFL(L+1))GO TO 273 ( 171) PBIN=PTIN ( 172) PHBI=PHTI ( 173) ENDDO ( 174) ! ( 175) 273 DPOSP=(PHTI-DFL(L+1))/TRTV ( 176) PRBIN=(1.+DPOSP)/(1.-DPOSP)*PTIN ( 177) ( 178) !*** ( 179) !*** FIND P AT THE REFERENCE INTERFACE GEOPOTENTIAL AT THE TOP ( 180) !*** ( 181) PBIN=PT+PD(I,J,n) ( 182) PHBI=DFL(LMST+1) ( 183) ! ( 184) DO LI=LMST,1,-1 ( 185) PTIN=PBIN-DETA(LI)*PD(I,J,n)*RES(I,J,n) ( 186) TRTV=2.*R*T(I,J,n,LI)*(1.+0.608*Q(I,J,n,LI)) ( 187) PHTI=PHBI+TRTV*(PBIN-PTIN)/(PBIN+PTIN) ( 188) IF(PHTI.GE.DFL(L))GO TO 275 ( 189) PBIN=PTIN ( 190) PHBI=PHTI ( 191) ENDDO ( 192) ! ( 193) 275 DPOSP=(PHTI-DFL(L))/TRTV ( 194) PRTIN=(1.+DPOSP)/(1.-DPOSP)*PTIN ( 195) ! ( 196) TTV(I,J,n)=(DFL(L)-DFL(L+1))/(2.*R)*(PRBIN+PRTIN)/(PRBIN-PRTIN) ( 197) ENDIF ( 198) ENDDO ( 199) ENDDO ( 200) ENDDO ( 201) ( 202) !---------------------------------------------------------------- ( 203) KMM=KMNTM ( 204) !---------------------------------------------------------------- ( 205) !*** ( 206) !*** HERE IS THE RELAXATION LOOP ( 207) !*** ( 208) !---------------------------------------------------------------- ( 209) DO 285 N1=1,NRLX ( 210) ! ( 211) CALL bocoh(TTV,im,jm,nm) ! Exchange haloes ( 212) ! ( 213) DO 280 KM=1,KMM ( 214) I=IMNT(KM) ( 215) J=JMNT(KM) PGF90 (Version 12.8) 08/30/2020 00:27:52 page 7 ( 216) n=nMNT(KM) ( 217) if(isCorner2(i,j,n,im,jm,nm).eq.0)then ( 218) TTV(I,J,n)=AD05*(4.*(TTV(I,J-1,n)+TTV(I+1,J,n) & ( 219) +TTV(I,J+1,n)+TTV(I-1,J,n)) & ( 220) +TTV(I-1,J-1,n) +TTV(I-1,J+1,n) & ( 221) +TTV(I+1,J-1,n) +TTV(I+1,J+1,n)) & ( 222) -CFT0*TTV(I,J,n) ( 223) else ( 224) TTV(i,j,n)=0.25*(TTV(I,J-1,n)+TTV(I+1,J,n) & ( 225) +TTV(I,J+1,n)+TTV(I-1,J,n)) ( 226) endif ( 227) 280 CONTINUE ( 228) ! ( 229) 285 CONTINUE ( 230) !---------------------------------------------------------------- ( 231) ! ( 232) if(l.eq.36)then ( 233) ! print *,'t(30,101,3),',t(30,101,3,36),ttv(30,101,3) ( 234) endif ( 235) ( 236) DO 290 KM=1,KMM ( 237) I=IMNT(KM) ( 238) J=JMNT(KM) ( 239) n=nMNT(KM) ( 240) if(l.eq.36)then ( 241) ! print *,i,j,n,km ( 242) endif ( 243) T(I,J,n,L)=TTV(I,J,n) ( 244) 290 CONTINUE ( 245) ! ( 246) 300 CONTINUE ( 247) ! VALUES FOR IMNT AND JMNT ARE FOR LAYER LM - THIS IS WHAT WE WANT ( 248) KMM=KMNTM ( 249) ! ( 250) DO 320 KM=1,KMM ( 251) I=IMNT(KM) ( 252) J=JMNT(KM) ( 253) n=nMNT(KM) ( 254) LMAP1=LMH(I,J,n)+1 ( 255) PBIN=PT+PD(I,J,n) ( 256) ! ( 257) DO L=LMAP1,LM ( 258) PTIN=PBIN ( 259) DPOSP=(DFL(L)-DFL(L+1))/(2.*R*T(I,J,n,L)) ( 260) PBIN=(1.+DPOSP)/(1.-DPOSP)*PTIN ( 261) if((i.eq.35.or.i.eq.36).and.j.eq.58.and.n.eq.3)then ( 262) ! print *,i,j,n,l,pd(i,j,n),ptin,dposp ( 263) ! stop ( 264) endif ( 265) ( 266) ENDDO ( 267) ! ( 268) PSLP(I,J,n)=PBIN ( 269) if(i.eq.70.and.j.eq.186.and.n.eq.4)then ( 270) ! print *,i,j,n,lmap1,lm ( 271) ! stop ( 272) endif ( 273) if(pslp(i,j,n).gt.123000)then PGF90 (Version 12.8) 08/30/2020 00:27:52 page 8 ( 274) ! print *,i,j,n,pslp(i,j,n) ( 275) ! stop ( 276) endif ( 277) 320 CONTINUE ( 278) ! stop ( 279) ( 280) else ( 281) ( 282) write(6,*) 'doing standard reduction!!!' ( 283) DO 410 n=1,nM ( 284) DO 410 J=1,JM ( 285) DO 410 I=1,IM ( 286) IF(FIS(I,J,n).GE.1.)THEN ( 287) LMA=LMH(I,J,n) ( 288) ALPP1=ALOG(PDSL1(I,J,n)*ETA(LMA+1)+PT) ( 289) SLOP=0.0065*ROG*T(I,J,n,LMA) ( 290) IF(SLOP.LT.0.50)THEN ( 291) SLPP=ALPP1+FIS(I,J,n)/(R*T(I,J,n,LMA)) ( 292) ELSE ( 293) TTT=-(ALOG(PDSL1(I,J,n)*ETA(LMA)+PT)+ALPP1) & ( 294) *SLOP*0.50+T(I,J,n,LMA) ( 295) SLPP=(-TTT+SQRT(TTT*TTT+2.*SLOP* & ( 296) (FIS(I,J,n)/R+ & ( 297) (TTT+0.50*SLOP*ALPP1)*ALPP1)))/SLOP ( 298) ENDIF ( 299) PSLP(I,J,n)=EXP(SLPP) ( 300) ENDIF ( 301) 410 CONTINUE ( 302) endif ( 303) ( 304) do 440 n=1,nm ( 305) DO 440 J=1,JM ( 306) DO 440 I=1,IM ( 307) SLPX(I,J,n)=PSLP(I,J,n) ( 308) 440 CONTINUE ( 309) !C ( 310) !C ( 311) !C!$OMP parallel do private(ihh2) ( 312) do 460 n=1,nm ( 313) DO 460 J=2,JM-1 ( 314) DO 460 I=2,im-1 ( 315) !C ( 316) !C*** EXTRA AVERAGING UNDER MOUNTAINS TAKEN OUT, FM, MARCH 96 ( 317) !C ( 318) SLPX(I,J,n)=0.125*(PSLP(I,J-1,n)+PSLP(I,J+1,n) & ( 319) +PSLP(I+1,J,n)+PSLP(I-1,J,n) & ( 320) +4.*PSLP(I,J,n)) ( 321) 460 CONTINUE ( 322) !C ( 323) !C!$OMP parallel do ( 324) do n=1,nm ( 325) DO J=1,JM ( 326) DO I=1,IM ( 327) PSLP(I,J,n)=SLPX(I,J,n) ( 328) ENDDO ( 329) ENDDO ( 330) ENDDO ( 331) PGF90 (Version 12.8) 08/30/2020 00:27:52 page 9 ( 332) end subroutine slp