subroutine eta2p USE mod_vrbls USE mod_extra USE mod_masks USE mod_pvrbls USE mod_graph implicit none include 'param_o.h' !!JNT include 'extra.comm' include 'mapot.comm' include 'const.h' include 'dynam.comm' include 'params' integer,parameter::im_jm_nm=(im+2)*(jm+2)*nm real,parameter::g=9.80616,RD1=287.04,GAMMA=6.5E-3, & RGAMOG=RD1*GAMMA/G,PQ0=379.90516,A2=17.2693882 & , A3=273.16,A4=35.86 integer::nhold,lp,i,j,n,lsl,nn,l,lmb,il !!JNT real,dimension(0:im+1,0:jm+1,nm)::tsl,qsl,fsl,q2sl !!JNT real,dimension(0:im+1,0:jm+1,nm,lm)::iw !!JNT integer,dimension(0:im+1,0:jm+1,nm)::nl1x !!JNT integer,dimension(im_jm_nm)::ihold,jhold,nnhold real, allocatable, dimension(:,:,:)::tsl,qsl,fsl,q2sl real, allocatable, dimension(:,:,:,:)::iw integer, allocatable, dimension(:,:,:)::nl1x integer, allocatable, dimension(:)::ihold,jhold,nnhold real::pnl1,B,fac,pu,tu,tabv,pl,tl,tabo,ahf,tblo,petau,alpetu, & petal,alpetl,alpet2,fact,alpet1,trf,zu,qu,ai,bi,tmt0,tmt15, & qw,qi,qint,qsat,utim,climit,lml,hh,tkl,qkl,cwmkl,pp, & u00kl,fiq,iwu,qabv,rhu,bq,ahfq,ql,iwl,rhl,qblo,q2a,bq2, & q2b,ahfq2 !!JNT real,dimension(0:im+1,0:jm+1,nm)::alpetux,alpet2x,usl,vsl real, allocatable, dimension(:,:,:)::alpetux,alpet2x,usl,vsl ! print *,'u,v,',u(10,10,6,6),v(10,10,6,6) ! print *,'us,vs,',u(10,10,6,6)*u2us(10,10,6)+v2us(10,10,6)*v(10,10,6,6) ! print *,'us,vs,',u(10,10,6,6)*u2vs(10,10,6)+v2vs(10,10,6)*v(10,10,6,6) UTIM=1. CLIMIT =1.0E-20 ALLOCATE(tsl(0:im+1,0:jm+1,nm)) ALLOCATE(qsl(0:im+1,0:jm+1,nm)) ALLOCATE(fsl(0:im+1,0:jm+1,nm)) ALLOCATE(q2sl(0:im+1,0:jm+1,nm)) ALLOCATE(iw(0:im+1,0:jm+1,nm,lm)) ALLOCATE(nl1x(0:im+1,0:jm+1,nm)) ALLOCATE(ihold(im_jm_nm)) ALLOCATE(jhold(im_jm_nm)) ALLOCATE(nnhold(im_jm_nm)) ALLOCATE(alpetux(0:im+1,0:jm+1,nm)) ALLOCATE(alpet2x(0:im+1,0:jm+1,nm)) ALLOCATE(usl(0:im+1,0:jm+1,nm)) ALLOCATE(vsl(0:im+1,0:jm+1,nm)) do n=1,nm DO J=1,jm DO I=1,IM IW(I,J,n,1)=0. ENDDO ENDDO ENDDO !C DO L=2,LM do n=1,nm DO J=1,jm DO I=1,IM LML=LM-LMH(I,J,n) HH=HTM(I,J,n,L) TKL=T(I,J,n,L) QKL=Q(I,J,n,L) CWMKL=W(I,J,n,L) TMT0=(TKL-273.16)*HH TMT15=AMIN1(TMT0,-15.)*HH PP=PDSL(I,J,n)*AETA(L)+PT QW=HH*PQ0/PP*EXP(HH*A2*(TKL-A3)/(TKL-A4)) QI=QW*(1.+0.01*AMIN1(TMT0,0.)) U00KL=U00(I,J,n)+UL(L+LML)*(0.95-U00(I,J,n))*UTIM !C IF(TMT0.LT.-15.0)THEN FIQ=QKL-U00KL*QI IF(FIQ.GT.0..OR.CWMKL.GT.CLIMIT) THEN IW(I,J,n,L)=1. ELSE IW(I,J,n,L)=0. ENDIF !C IF(TMT0.GE.0.0)IW(I,J,n,L)=0. IF(TMT0.LT.0.0.AND.TMT0.GE.-15.0)THEN IW(I,J,n,L)=0. IF(IW(I,J,n,L-1).EQ.1.0.AND.CWMKL.GT.CLIMIT)IW(I,J,n,L)=1. ENDIF endif enddo enddo enddo enddo open(1,file='tmpfld.dat',form='unformatted') lsl=lsm DO 310 LP=1,LSL NHOLD=0 TRF=2.*ALSL(LP) do 125 n=1,nm DO 125 J=1,jm DO 125 I=1,IM ! TSL(I,J,n)=-1.E6 QSL(I,J,n)=-1.E6 FSL(I,J,n)=-1.E6 ! !*** LOCATE VERTICAL INDEX OF MODEL INTERFACE JUST BELOW !*** THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING. ! DO 115 L=2,LM ! if(i.eq.25.and.j.eq.22.and.n.eq.3)then ! print *,l,alpint(i,j,n,l),alsl(lp),alpint(i,j,n,lp1) ! endif IF(ALPINT(I,J,n,L).GE.ALSL(LP))THEN NL1X(I,J,n)=L NHOLD=NHOLD+1 IHOLD(NHOLD)=I JHOLD(NHOLD)=J nnhold(nhold)=n GO TO 125 ENDIF 115 CONTINUE nl1x(i,j,n)=l NHOLD=NHOLD+1 IHOLD(NHOLD)=I JHOLD(NHOLD)=J nnhold(nhold)=n ! 125 CONTINUE DO 220 NN=1,NHOLD I=IHOLD(NN) J=JHOLD(NN) n=nnHOLD(NN) PNL1=PINT(I,J,n,NL1X(I,J,n)) IF(NL1X(I,J,n).EQ.1)THEN !--------------------------------------------------------------------- !*** EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL !--------------------------------------------------------------------- ! PU=PINT(I,J,n,2) ZU=ZINT(I,J,n,2) TU=0.5*(T(I,J,n,1)+T(I,J,n,2)) QU=0.50*(Q(I,J,n,1)+Q(I,J,n,2)) IWU=0.5*(IW(I,J,n,1)+IW(I,J,n,2)) TABV=TU*(SPL(LP)/PU)**RGAMOG TMT0=TU-273.16 TMT15=AMIN1(TMT0,-15.) AI=0.008855 BI=1. IF(TMT0.LT.-20.)THEN AI=0.007225 BI=0.9674 ENDIF QW=PQ0/PU*EXP(A2*(TU-A3)/(TU-A4)) QI=QW*(BI+AI*AMIN1(TMT0,0.)) QINT=QW*(1.-0.00032*TMT15*(TMT15+15.)) IF(TMT0.LT.-15.)THEN QSAT=QI ELSEIF(TMT0.GE.0.)THEN QSAT=QINT ELSE IF(IWU.GT.0.0) THEN QSAT=QI ELSE QSAT=QINT ENDIF ENDIF !CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT !C DELETE THIS LINE TO SWITCH BACK TO RH VS ICE QSAT=QW !CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT RHU =QU/QSAT !C IF(RHU.GT.1)THEN RHU=1 QU =RHU*QSAT ENDIF IF(RHU.LT.0.01)THEN RHU=0.01 QU =RHU*QSAT ENDIF TMT0=TABV-273.16 TMT15=AMIN1(TMT0,-15.) AI=0.008855 BI=1. IF(TMT0.LT.-20.)THEN AI=0.007225 BI=0.9674 ENDIF QW=PQ0/SPL(LP)*EXP(A2*(TABV-A3)/(TABV-A4)) QI=QW*(BI+AI*AMIN1(TMT0,0.)) QINT=QW*(1.-0.00032*TMT15*(TMT15+15.)) IF(TMT0.LT.-15.)THEN QSAT=QI ELSEIF(TMT0.GE.0.)THEN QSAT=QINT ELSE IF(IWU.GT.0.0) THEN QSAT=QI ELSE QSAT=QINT ENDIF ENDIF !CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT !C DELETE THIS LINE TO SWITCH BACK TO RH VS ICE QSAT=QW !CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT QABV =RHU*QSAT QABV =AMAX1(1.e-12,QABV) B =TABV BQ =QABV FAC =0. AHF =0. AHFQ =0. Q2A =0.50*(Q2(I,J,n,1)+Q2(I,J,n,2)) BQ2 =Q2A ELSEIF(NL1X(I,J,n).EQ.LP1)THEN !--------------------------------------------------------------------- !*** EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND) !--------------------------------------------------------------------- ! PL=PINT(I,J,n,LM-1) TL=0.5*(T(I,J,n,LM-2)+T(I,J,n,LM-1)) QL=0.5*(Q(I,J,n,LM-2)+Q(I,J,n,LM-1)) IWL=0.50*(IW(I,J,n,LM-2)+IW(I,J,n,LM-1)) TBLO=TL*(SPL(LP)/PL)**RGAMOG TMT0=TL-273.16 TMT15=AMIN1(TMT0,-15.) AI=0.008855 BI=1. IF(TMT0.LT.-20.)THEN AI=0.007225 BI=0.9674 ENDIF QW=PQ0/PL*EXP(A2*(TL-A3)/(TL-A4)) QI=QW*(BI+AI*AMIN1(TMT0,0.)) QINT=QW*(1.-0.00032*TMT15*(TMT15+15.)) IF(TMT0.LT.-15.)THEN QSAT=QI ELSEIF(TMT0.GE.0.)THEN QSAT=QINT ELSE IF(IWL.GT.0.0) THEN QSAT=QI ELSE QSAT=QINT ENDIF ENDIF !CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT !C DELETE THIS LINE TO SWITCH BACK TO RH VS ICE QSAT=QW !CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT RHL=QL/QSAT IF(RHL.GT.1)THEN RHL=1 QL =RHL*QSAT ENDIF IF(RHL.LT..01)THEN RHL=.01 QL =RHL*QSAT ENDIF TMT0=TBLO-273.16 TMT15=AMIN1(TMT0,-15.) AI=0.008855 BI=1. IF(TMT0.LT.-20.)THEN AI=0.007225 BI=0.9674 ENDIF QW=PQ0/SPL(L) & *EXP(A2*(TBLO-A3)/(TBLO-A4)) QI=QW*(BI+AI*AMIN1(TMT0,0.)) QINT=QW*(1.-0.00032*TMT15*(TMT15+15.)) IF(TMT0.LT.-15.)THEN QSAT=QI ELSEIF(TMT0.GE.0.)THEN QSAT=QINT ELSE IF(IWL.GT.0.0) THEN QSAT=QI ELSE QSAT=QINT ENDIF ENDIF !CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT !C DELETE THIS LINE TO SWITCH BACK TO RH VS ICE QSAT=QW !CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT QBLO =RHL*QSAT QBLO =AMAX1(1.e-12,QBLO) B =TBLO BQ =QBLO FAC =0. AHF =0. AHFQ =0. Q2A =0.50*(Q2(I,J,n,LMH(I,J,n)-1)+Q2(I,J,n,LMH(I,J,n))) BQ2 =Q2A ELSE !--------------------------------------------------------------------- !*** INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS !--------------------------------------------------------------------- ! B =T(I,J,n,NL1X(I,J,n)) BQ =Q(I,J,n,NL1X(I,J,n)) FAC =2.*ALOG(PT+PDSL(I,J,n)*AETA(NL1X(I,J,n))) AHF =(B-T(I,J,n,NL1X(I,J,n)-1))/ & (ALPINT(I,J,n,NL1X(I,J,n)+1)-ALPINT(I,J,n,NL1X(I,J,n)-1)) AHFQ =(BQ-Q(I,J,n,NL1X(I,J,n)-1))/ & (ALPINT(I,J,n,NL1X(I,J,n)+1)-ALPINT(I,J,n,NL1X(I,J,n)-1)) Q2B =0.50*(Q2(I,J,n,NL1X(I,J,n)-1)+Q2(I,J,n,NL1X(I,J,n))) !C IF(NL1X(I,J,n).GT.2)THEN Q2A=0.50*(Q2(I,J,n,NL1X(I,J,n)-2)+Q2(I,J,n,NL1X(I,J,n)-1)) ELSE Q2A=Q2B ENDIF !C BQ2=Q2B*HTM(I,J,n,NL1X(I,J,n)) AHFQ2=(BQ2-Q2A*HTM(I,J,n,NL1X(I,J,n)-1))/ & (ALPINT(I,J,n,NL1X(I,J,n)+1)-ALPINT(I,J,n,NL1X(I,J,n)-1)) if(i.eq.25.and.j.eq.22.and.n.eq.3)then ! print *,'tsl,',lp,nl1x(i,j,n),ahf,ALPINT(I,J,n,NL1X(I,J,n)+1) & ! ,ALPINT(I,J,n,NL1X(I,J,n)-1) endif ENDIF TSL(I,J,n)=B+AHF*(TRF-FAC) QSL(I,J,n)=BQ+AHFQ*(TRF-FAC) QSL(I,J,n)=AMAX1(QSL(I,J,n),1.e-12) Q2SL(I,J,n)=BQ2+AHFQ2*(TRF-FAC) Q2SL(I,J,n)=AMAX1(Q2SL(I,J,n),0.) FSL(I,J,n)=(PNL1-SPL(LP))/(SPL(LP)+PNL1) & *((ALSL(LP)+ALPINT(I,J,n,NL1X(I,J,n))-FAC)*AHF+B)*Rd*2. & +ZINT(I,J,n,NL1X(I,J,n))*G if(i.eq.25.and.j.eq.22.and.n.eq.3)then ! print *,'tsl,',lp,tsl(i,j,n),b,ahf,trf,fac,nl1x(i,j,n) endif 220 CONTINUE do 281 n=1,nm DO 281 J=1,jm DO 281 I=1,IM !NOTE !NOTE 29 JANUARY 1993, RUSS TREADON. !NOTE - AS FOR THE OTHER FIELDS WE INTERPOLATE ONLY !NOTE BETWEEN THE FAL AND THE MODEL TOP. BELOW !NOTE SURFACE VALUES ARE FAL VALUES. ! LMB = LMV(I,J,n) ! PETAU=PT+PDVP1(I,J,n)*ETA(1) ALPETU=ALOG(PETAU) DO 280 IL=2,LMB PETAL=PT+PDVP1(I,J,n)*ETA(IL) ! PETAU=PT+PDVP1(I,J)*ETA(IL-1) ALPETL=ALOG(PETAL) ! ALPETU=ALOG(PETAU) ALPET2=SQRT(0.5E0*(ALPETL*ALPETL+ALPETU*ALPETU)) ! ! SEARCH FOR HIGHEST MID-LAYER ETA SURFACE (NOT SUBMERGED) ! THAT IS BELOW THE GIVEN STANDARD PRESSURE LEVEL. IF(ALSL(Lp).LT.ALPET2)THEN NL1X(I,J,n)=IL-1 ALPETUX(I,J,n)=ALPETU ALPET2X(I,J,n)=ALPET2 GO TO 281 ENDIF ! If we arent on the last iterate of the 280 loop, reset PETAU and ALPETU if ( il .eq. lmb ) goto 280 PETAU=PETAL ALPETU=ALPETL 280 CONTINUE NL1X(I,J,n)=LMB+1 ALPETUX(I,J,n)=ALPETU ALPET2X(I,J,n)=ALPET2 281 CONTINUE ! ! BELOW GROUND USE FAL WINDS. ! !$omp parallel do !$omp private(alpet1,alpetl,alpetu,fact,petau) do 290 n=1,nm DO 290 J=1,jm DO 290 I=1,IM IF(NL1X(I,J,n).GT.LMV(I,J,n))THEN USL(I,J,n)=U(I,J,n,LMV(I,J,n)) VSL(I,J,n)=V(I,J,n,LMV(I,J,n)) ! ! IF REQUESTED PRESSURE LEVEL IS NOT BELOW THE LOCAL GROUND ! THEN WE HAVE TWO POSSIBILITIES. IF THE REQUESTED PRESSURE ! LEVEL IS BETWEEN THE LOCAL SURFACE PRESSURE AND TOP OF ! MODEL PRESSURE, VERTICALLY INTERPOLATE BETWEEN NEAREST ! BOUNDING ETA LEVELS TO GET THE WIND COMPONENTS. IF THE ! REQUESTED PRESSURE LEVEL IS ABOVE THE MODEL TOP, USE ! CONSTANT EXTRAPOLATION OF TOP ETA LAYER (L=1) WINDS. ! ELSE IF(NL1X(I,J,n).GT.1)THEN ALPETL=ALPETUX(I,J,n) PETAU=PT+PDVP1(I,J,n)*ETA(NL1X(I,J,n)-1) ALPETU=ALOG(PETAU) ALPET1=SQRT(0.5*(ALPETL*ALPETL+ALPETU*ALPETU)) FACT=(ALPET2X(I,J,n)-ALSL(Lp))/(ALPET2X(I,J,n)-ALPET1) USL(I,J,n)=U(I,J,n,NL1X(I,J,n)) & +(U(I,J,n,NL1X(I,J,n)-1)-U(I,J,n,NL1X(I,J,n)))*FACT VSL(I,J,n)=V(I,J,n,NL1X(I,J,n))+(V(I,J,n,NL1X(I,J,n)-1) & -V(I,J,n,NL1X(I,J,n)))*FACT ELSE USL(I,J,n)=U(I,J,n,NL1X(I,J,n)) VSL(I,J,n)=V(I,J,n,NL1X(I,J,n)) ENDIF ! ! ALPET2 IS MID-LAYER ETA SURFACE JUST BELOW STANDARD PRESSURE ! LEVEL AND ALPET1 IS DASHED ETA SURFACE JUST ABOVE. ! NOTE THAT IF THE STANDARD PRESSURE SURFACE IS SUBMERGED, THEN ! ALPET2 AND ALPET1 ARE THE LOWEST AND 2ND LOWEST MID-LAYER ! ETA SURFACES ABOVE THE TOPOGRAPHY (WITH OLDRD=.TRUE., ZJ). ! ENDIF ! print *,i,j,n,usl(i,j,n),vsl(i,j,n) 290 CONTINUE do j=1,jm do i=1,im ! print *,i,j,q2sl(i,j,3) enddo enddo ! print *,'usl,99,58,6,',usl(99,58,6),usl(99,59,6) ! print *,'vsl,99,58,6,',lp,vsl(99,58,6),vsl(99,59,6) call bocoh(fsl,im,jm,nm) call bocov(usl,vsl,im,jm,nm) ! print *,'usl,99,58,6,',usl(99,58,6),usl(99,59,6) ! print *,'vsl,99,58,6,',lp,vsl(99,58,6),vsl(99,59,6) write(1)fsl write(1)usl write(1)vsl write(1)qsl !DRAGAN 23.09.'10 write(1)tsl !!!avg2017 write(1)q2sl !!! write(1)w !!! write(1)dwdt !non-hydrostatic 2016 !!! write(1)z !! write(1)tshltr !! write(1)qshltr ! write(1)plm 310 continue close(1) DEALLOCATE(tsl,qsl,fsl,q2sl,iw) DEALLOCATE(nl1x,ihold,jhold,nnhold) DEALLOCATE(alpetux,alpet2x,usl,vsl) end subroutine eta2p