subroutine bulkadjust(iee) !----------------------------------------------------------------------- !----------------------------------------------------------------------- !16 January 2002 Tuleya, Khain and Pokrovsky ! ADJUST handles all of the adjustments due to condensation. the ! computations are not done in the hole. the final adjusted ! temperature and moisture are put into "tf" and "rf". the ! actual tendencys produced are stored for diagnostic purposes ! "tdt3" and tdr3". !----------------------------------------------------------------------- ! use allocate_mod use module_mp_lin include 'RESOLUTION.h' include 'PARAMETERS.h' include 'STDUNITS.h' include 'NUM.h' include 'AREAS.h' include 'BKINFO.h' include 'COMMUNICATE.h' include 'CONSLEV.h' include 'CONMLEV.h' include 'CLCLE.h' include 'ESTAB.h' include 'GDINFO.h' include 'FILEIF.h' include 'FILEPC.h' include 'FLAGS.h' include 'FLEXK.h' include 'LIMIT.h' include 'BINI1I2.h' ! include 'BULKCHUNK.h' include 'QLOGS.h' include 'RADDAT.h' include 'SET.h' include 'STORE.h' include 'STORJ.h' include 'TIME.h' include 'TEND.h' include 'CONVEC.h' include 'NUNIT.h' common /detrain/ sasm(kmax,imx,2) common /idiag/ idiagn common /retain/ fac1(imx,jmx,kmax,3), * fac2(70,70,kmax,3), * fac3(94,94,kmax,3) common /isch/ ischeme,imom common /micros/ ibulk common /tjunk/ uold(kmax,imx),vold(kmax,imx), * told(kmax,imx),rold(kmax,imx) cc * ta(kmax,imx),ra(kmax,imx) !----------------------------------------------------------------------- ! internal variables !----------------------------------------------------------------------- real, dimension(imx) :: precm real, dimension(imx) :: precc real, dimension(kmax,imx) :: tcon real, dimension(kmax,imx) :: rcon real, dimension(kmax,imx) :: udif real, dimension(kmax,imx) :: vdif real, dimension(kmax,imx) :: tdif real, dimension(kmax,imx) :: rdif real, dimension(kmax,imx) :: tmew real, dimension(kmax,imx) :: rmew real, dimension(kmax,i1:i2) :: ra real, dimension(kmax,i1:i2) :: ta real, dimension(imx) :: ttlprc real, dimension(i1:i2,kmax+1) :: zhalf real, dimension(i1:i2,kmax) :: zfull !----------------------------------------------------------------------- ! BULK variables !----------------------------------------------------------------------- real, dimension(iee ,kmax) :: ta_hum real, dimension(iee ) :: slmsk real, dimension(iee ) :: aprec real, dimension(iee ) :: sr real, dimension(iee ,kmax) :: qvr real, dimension(iee ,kmax) :: qcl real, dimension(iee ,kmax) :: qrn real, dimension(iee ,kmax) :: qci real, dimension(iee ,kmax) :: F_Rime real, dimension(iee ,kmax) :: F_Ice real, dimension(iee ,kmax) :: F_Rain real, dimension(iee ,kmax) :: qcc real, dimension(iee ,kmax) :: ss real, dimension(iee ,kmax) :: tv real, dimension(iee ,kmax) :: tha real, dimension(iee ,kmax) :: pi real, dimension(iee ,kmax) :: rho real, dimension(iee ,kmax) :: dz8w real, dimension(iee ,kmax) :: th real, dimension(iee ,kmax) :: g_cond real, dimension(iee ,kmax) :: g_evap real, dimension(iee ,kmax) :: g_dep real, dimension(iee ,kmax) :: g_subl real, dimension(iee ,kmax) :: g_freeze real, dimension(iee ,kmax) :: g_melt ! real, dimension(iee) :: pss real, dimension(iee ,kmax) :: z real, dimension(iee ,kmax) :: psq real, dimension(iee ,kmax) :: dqs real, dimension(iee ,kmax) :: rhc real, dimension(iee ) :: RAINNC real, dimension(iee) :: ht real, dimension(iee) :: xncw real, dimension(imx) :: rain logical lprnt ! INTEGER ids,ide, jds,jde, kds,kde , * ims,ime, jms,jme, kms,kme , * its,ite, jts,jte, kts,kte , * P_QI,P_QS,P_QG, * I_FIRST_SCALAR ! character*10 routine ! !----------------------------------------------------------------------- ! code !----------------------------------------------------------------------- epsln = 1.0e-20 tbat = 3.0e+02 arcppp=2./7. c c N1 determine the droplet concentration (n/cm3) c Which we have set to 100 for maritime air c tmnus = -99.9 c if(ttime(ngd) .EQ. 0.0) tbat = deltat j = IFIX(tjlof(i1)) if(iacp .EQ. 0 .AND. l .EQ. 1) then if(j .EQ. 1 .AND. ngd .EQ. 1) write(stdout,10) nstep 10 format(' PRECIPITATION RESET AT TIME STEP (NSTEP) = ', i5) do i = i1,i2 acp6f(i) = 0.0 enddo endif !----------------------------------------------------------------------- ! initialize arrays !----------------------------------------------------------------------- do i = i1,i2 do k = 1,kmax told(k,i) = tf(k,i) rold(k,i) = rf(k,i,1) uold(k,i) = uf(k,i) vold(k,i) = vf(k,i) ta (k,i) = tf(k,i) ra (k,i) = rf(k,i,1) enddo enddo !----------------------------------------------------------------------- ! check that temperature range does not exceed estab range !----------------------------------------------------------------------- do i = i1,i2 do k = 1,kmax if(ta(k,i) .GT. 375.1) go to 30 enddo enddo do i = i1,i2 do k = 1,kmax if(ta(k,i) .LT. 100.2) go to 30 enddo enddo go to 50 30 continue print*,'temperatures out of range' write(6,35) k, i 35 format(' TEMPERATURE POINT IN "BULKADJ" REQUIRE TERMINATION STATUS & ', 2(I5,2X)) write(6,40) j, ngd, (ta(k,i),k=1,kmax) 40 format(' TEMPERATURES NOW IN ROW ', I5,' GRID NUMBER ', I5, /, &(4X,10E11.4)) c thour = ttime(nnest)/3600.0 ihours = IFIX(thour) ihours = ihours - 1 print*,' proceed with postprocessing' rewind kstats_6 do iii = 1 ,ihours, 6 print *,'this is iii :',iii,ihours read (kstats_6,390) hour, poslon, poslat, psmin, wsfmk write(6,390) hour, poslon, poslat, psmin, wsfmk write (kstats_6,390) hour, poslon, poslat, psmin, wsfmk 390 format('HOUR:',F5.1,' LONG:',F8.2,' LAT:',F7.2, & ' MIN PRESS (hPa): ',F7.2,' MAX SURF WIND (KNOTS):',F6.2) enddo print*,'writing final stats marker at hour: ',ihours write(kstats_6,65) tmnus, tmnus, tmnus write(kstat1,*) ihours routine = 'bulkadjust' call MPI_CLOSE(1,routine) 65 format(5X, F5.1, 7X, F8.2, 6X, F7.2) 50 continue !----------------------------------------------------------------------- ! only do convective adjustment during corrector step !----------------------------------------------------------------------- if(dift(1,ngd).GE. 23 .OR. dift(2,ngd) .GE. 23) then toff = 1.0 else toff = 0.0 endif if(toff .EQ. 1 .AND. (I1.EQ.1.OR.I2.EQ.imax))go to 420 if(toff .EQ. 1 .AND. (j.LE. 2+ngr .OR. J .GE.jmax-1-ngr)) * go to 420 ! print *,'itr in bulkajustBEGIN',itr if(itr .EQ. 0 .AND. l .EQ. 0) go to 420 ! ! now do other water variables. do in relative i space . ! check if chunk size is too big ! ! i2mi1p = i2 -i1 +1 ! if(i2mi1p .GT. nx) then ! write(stdout,45) ngd, i1, i2, nx !45 format(' MODEL CHUNK > MICRO CHUNK FOR GRID ',3i4,' i1,i2,nx') ! routine = 'bulkadjust' ! call MPI_CLOSE(1,routine) ! endif c c do SAS convective paramaterization after micro-physics c do k = 1,kmax do i = i1,i2 tmew(k,i) = ta(k,i) rmew(k,i) = ra(k,i) tcon(k,i) = ta(k,i) rcon(k,i) = ra(k,i) udif(k,i) = 0.0 vdif(k,i) = 0.0 tdif(k,i) = 0.0 rdif(k,i) = 0.0 enddo enddo c !----------------------------------------------------------------------- if(ischeme .EQ. 1) then call Setconv (uf,vf,tcon,rcon, psf, i1, i2, * inc,j ,tdif, rdif, udif, vdif, precc,imom) do i = i1,i2 do k = 1,kmax uf(k,i) = uold(k,i) + udif(k,i) vf(k,i) = vold(k,i) + vdif(k,i) ta(k,i) = tmew(k,i) + tdif(k,i) ra(k,i) = rmew(k,i) + rdif(k,i) enddo enddo else call setcone (ua, va, ta, ra, rs, pnew, i1, i2, * inc,j ,tdif, rdif ,udif,vdif, precc,imom) do i = i1,i2 do k = 1,kmax uf(k,i) = uold(k,i) + udif(k,i) * deltat vf(k,i) = vold(k,i) + vdif(k,i) * deltat ta(k,i) = tmew(k,i) + tdif(k,i) * deltat ra(k,i) = rmew(k,i) + rdif(k,i) * deltat enddo enddo endif c if(ngd .NE. nnest) then do i = i1,i2 if(j.GT. js .AND. j .LT. jn .AND. i .GT. iw .AND. i .LT. ie)then do k = 1,kmax uf(k,i) = uold(k,i) vf(k,i) = vold(k,i) ta(k,i) = told(k,i) ra(k,i) = rold(k,i) enddo endif enddo endif !----------------------------------------------------------------------- ! calculate variables needed for GODDARD-GFDL bulk package !-------------------------------------------------------------------- ! print *, 'start doing bulk' call height (ta, psf, zf, q, qmh, kmax,i1,i2,imx ,zhalf,zfull) ! ! ! constants used for Lin scheme ! n_moist = 6 I_FIRST_SCALAR = 0 P_QV = 1 P_QC = 2 P_QR = 3 P_QI = 4 P_QS = 5 P_QG = 6 poo = 1.0e5 rhoair0 = 1.28 rhowater = 1000. rhosnow = 100. EPS=0.62197 Rsd = 287.0 Rsv = 461.6 CPP = 1004.5 CV = 717.5 GG = .01 * GGG RCP = Rsd/CPP RCPM= -Rsd/CPP XLS = 2.85E6 XLV = 2.5E6 XLF = 3.50E5 SVP1=0.6112 SVP2=17.67 SVP3=29.65 SVPT0=273.15 EP_1=Rsv/Rsd-1. EP_2=Rsd/Rsv ! c remove the negative condensate c c determine the highest level you want to borrow from c istart=12 do isp = 2 , 4 c do i = i1 , i2 eps = 9.99999996E-13 sumn = 0.0 qsump = 0.0 tt = 0.0 c c determine the amount of condensate you need to borrow c do k = 1 , kmax if(rf(k,i,isp) .LT. 0.0) then tt = 1.0 sumn = sumn + (rf(k,i,isp)-eps) * dq(k) endif enddo c c borrow the condensate from levels that have enough positive amount c testa = -sumn*1.2 do k = istart , kmax if(rf(k,i,isp) .GT. testa)then qsump = qsump + dq(k) endif enddo do k = istart , kmax if(rf(k,i,isp) .GT. testa) then rf(k,i,isp) = rf(k,i,isp) + sumn*dq(k)/qsump endif enddo c c put the borrowed condensate into the levels that are negative c do k = 1 , kmax if(rf(k,i,isp) .LT. 0.0) then rf(k,i,isp) = eps endif enddo c c proceed with the next i point c enddo c enddo c do i = i1,i2 nii = i -i1 +1 ht(nii) = .01 * zf(i)/ggg pss(nii) = .1 * psf(i) do k = 1,kmax kz=kmax + 1 - k ta_hum(nii,kz) = ta(k,i) qvr(nii,kz) = ra (k,i) z(nii,kz) = .01 * zfull(i,k)/ggg ! ! mks units is a factor of .1 from cgs for pressure : ! psq(nii,kz) = .1 * psf(i)*qv(k,i) c c c input combined condensate to the driveer c c 3d change c c set comnbined condensate to 0 since no longer used c qcc(nii,kz)=0.0 c c input actual micro-physics particles into driver c contribution from SAS also needs to be added here c qrn(nii,kz)=rf(k,i,3) qcl(nii,kz)=rf(k,i,2) + sasm(k,i,2) qci(nii,kz)=rf(k,i,4) + sasm(k,i,1) c c fractional change no longer used but kept in here c to generalize code if needed in the future... c if(ngd .EQ. 1) then F_Rime(nii,kz)= fac1(i,j,k,1) F_Ice(nii,kz) = fac1(i,j,k,2) F_Rain(nii,kz)= fac1(i,j,k,3) endif if(ngd .EQ. 2) then F_Rime(nii,kz)= fac2(i,j,k,1) F_Ice(nii,kz) = fac2(i,j,k,2) F_Rain(nii,kz)= fac2(i,j,k,3) endif if(ngd .EQ. 3) then F_Rime(nii,kz)= fac3(i,j,k,1) F_Ice(nii,kz) = fac3(i,j,k,2) F_Rain(nii,kz)= fac3(i,j,k,3) endif enddo enddo C C DOING THE FERRIER SCHEME >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>. if(ibulk.eq.2) then ! call ferrier scheme.. ! define needed constants in GFS units (mks).. GRAV =9.8062E+0 HVAP =2.5000E+6 HFUS =3.3358E+5 HSUB=HVAP+HFUS ! lprnt = .true. lprnt = .false. ipr = 5 ! redefine pressures in terms of Kpa which are used in GFS ! store rime factor in snow field for ferrier scheme..keep>1. do i = i1,i2 nii = i -i1 +1 do k = 1,kmax kz=kmax + 1 - k psq(nii,kz) = 1.e-4 * psf(i)*qv(k,i) dqs(nii,kz) = qmh(k+1) - qmh(k) c c F_Rime(nii,kz) = max(1.0, F_RIME(nii,kz)) enddo enddo ! set number concentration to 50/cm3 for test......mb... do i = i1,i2 nii = i -i1 +1 xncw(nii)=60. enddo ! set critical large scale rh (set to operational gfdl values) if(ngd .EQ. 1) then do i = i1,i2 nii = i -i1 +1 do k = 1,33 kz=kmax + 1 - k rhc(nii,kz) = .985 enddo rhc(nii, kmax+1-34) = .99 do k = 35 , kmax kz=kmax + 1 - k rhc(nii,kz) = 1.0 enddo enddo else do i = i1,i2 nii = i -i1 +1 do k = 1,33 kz=kmax + 1 - k rhc(nii,kz) = .99 enddo rhc(nii, kmax+1-34) = .995 do k = 35 , kmax kz=kmax + 1 - k rhc(nii,kz) = 1.0 enddo enddo endif c CALL gferdrive(iee,iee,kmax,deltat,psq,pss,dqs, & ta_hum,qvr,qcc,slmsk, & qci,qcl,qrn,F_Rime,F_Ice,F_Rain, & aprec,sr,grav,hvap,hsub,cpp,rhc,xncw, & proc_num,lprnt,ipr) endif ! end of ferrier scheme.... ! !----------------------------------------------------------------------- ! put new variables back in forecast file ! do in relative i space !----------------------------------------------------------------------- kend = 1 do i =i1,i2 nii = i -i1 +1 precm(i) = 100. * aprec(nii) do k = kend , kmax kz=kmax+1-k ta(k,i)=ta_hum(nii,kz) ra(k,i) =qvr(nii,kz) c c all three species stored for advection rf(k,i,2)=qcl(nii,kz) rf(k,i,3)=qrn(nii,kz) rf(k,i,4)=qci(nii,kz) c c output cloud water, rain water and ice for diagnostics at history write c if(ihsflg .EQ. 1 .AND. idiagn .EQ. 1) then if(ngd .EQ. 1) then stor1(k,i,j,6) = qcl(nii,kz) stor1(k,i,j,7) = qrn(nii,kz) stor1(k,i,j,8) = qci(nii,kz) endif if(ngd .EQ. 2) then stor2(k,i,j,6) = qcl(nii,kz) stor2(k,i,j,7) = qrn(nii,kz) stor2(k,i,j,8) = qci(nii,kz) endif if(ngd .EQ. 3) then stor3(k,i,j,6) = qcl(nii,kz) stor3(k,i,j,7) = qrn(nii,kz) stor3(k,i,j,8) = qci(nii,kz) endif endif c c store the factional amounts in memory: if(ngd .EQ. 1) then fac1(i,j,k,1)=F_Rime(nii,kz) fac1(i,j,k,2)=F_Ice(nii,kz) fac1(i,j,k,3)=F_Rain(nii,kz) endif if(ngd .EQ. 2) then fac2(i,j,k,1)=F_Rime(nii,kz) fac2(i,j,k,2)=F_Ice(nii,kz) fac2(i,j,k,3)=F_Rain(nii,kz) endif if(ngd .EQ. 3) then fac3(i,j,k,1)=F_Rime(nii,kz) fac3(i,j,k,2)=F_Ice(nii,kz) fac3(i,j,k,3)=F_Rain(nii,kz) endif enddo 999 continue enddo c if(ngd .LE. 2) then do i = i1,i2 if(j.GT. js .AND. j .LT. jn .AND. i .GT. iw .AND. i .LT. ie)then do k = 1,kmax uf(k,i) = uold(k,i) vf(k,i) = vold(k,i) ta(k,i) = told(k,i) ra(k,i) = rold(k,i) rf(k,i,2) = 0.0 uold(k,i) = uf(k,i) vold(k,i) = vf(k,i) enddo endif enddo endif ! check that temperature range does not exceed estab range !----------------------------------------------------------------------- do i = i1,i2 do k = 1,kmax if(ta(k,i) .GT. 375.) go to 930 enddo enddo do i = i1,i2 do k = 1,kmax if(ta(k,i) .LT. 100.2) go to 930 enddo enddo go to 950 930 continue print*,'bad temperatures' write(6,935) k, i 935 format(' TEMPERATURE POINT IN "BULKADJ" AT POINT (K,I) OUT OF RANGE & ', 2(I5,2X)) write(6,940) j, ngd, (ta(k,i),k=1,kmax) 940 format(' TEMPERATURES BAD IN ROW ', I5,' GRID NUMBER ', I5, /, &(4X,10E11.4)) c c thour = ttime(nnest)/3600.0 ihours = IFIX(thour) ihours = ihours - 1 print*,' proceed with postprocessing' rewind kstats_6 do iii = 1 ,ihours, 6 print *,'this is iii; ',iii,ihours read (kstats_6,390) hour, poslon, poslat, psmin, wsfmk write(6,390) hour, poslon, poslat, psmin, wsfmk write (kstats_6,390) hour, poslon, poslat, psmin, wsfmk enddo print*,'writing final stats marker at hour: ',ihours write(kstats_6,65) tmnus, tmnus, tmnus write(kstat1,*) ihours routine = 'bulkadjust' call MPI_CLOSE(1,routine) 950 continue !----------------------------------------------------------------------- ! estimate precipitation during deltat(delt) !----------------------------------------------------------------------- if(l .EQ. 0) go to 420 do i = i1,i2 ttlprc(i) = 0.0 enddo c c c shallow and deep convective precip stored in : precc(i) c c with detrained micro-phyics convective precip no longer moisture c difference c c do k = 1,kmax c do i = i1,i2 c ttlprc(i) = ttlprc(i) - rdif(k,i)*dqv(k,i) c enddo c enddo c factr = rdeltx*deltat/ggg c c c Add in convective and large-scale precipitation c do i = i1,i2 ttlprc(i) = precc(i) + precm(i) c acp6f (i) = acp6f(i) + ttlprc(i) cbug acprf (i) = acp6f(i) acprf(i) = acprf(i) + ttlprc(i) enddo if(ngd .EQ. 333 .AND. j .EQ. 30) then cc write(6,3444)(acp6f(i), i = i1 , i2) cc write(6,3445)(precl(i), i = i1 , i2) cc write(6,3446)(acprf(i), i = i1 , i2) 3444 format('sas ',(/,12e12.6) ) 3445 format('brad',(/,12e12.6) ) 3446 format('combined ',(/,12e12.6) ) endif if(mxmflg .NE. 0) then do i = i1,i2,inc if(j .GT. js .AND. j .LT. jn .AND. & i .GT. iw .AND. i .LT. ie) ttlprc(i) = 0. if(j .GT. js .AND. j .LT. jn .AND. & i .GT. iw .AND. i .LT. ie) then else ttlp(i,j) = ttlprc(i)/dstp endif enddo endif if(intflg .NE. 0) then do i = i1,i2,inc if(j .GT. js .AND. j .LT. jn .AND. & i .GT. iw .AND. i .LT. ie) then else areams(i,j,6,ngd) = areams(i,j,6,ngd) + ttlprc(i)*dlsf(i)* & deltat/dstp endif enddo endif do i = i1,i2 ra(1,i) = 1.0e-10 enddo ! ! adjust total condensate ............................ ! ! !----------------------------------------------------------------------- ! negative mixing ratio adjustments !----------------------------------------------------------------------- 420 continue do i = i1,i2 do k = 1,kmax - 2 if(ra(k,i) .LE. epsln) then ra(k+1,i) = ra(k+1,i) - ddqk(i,k)*(epsln - ra(k,i)) ra(k+2,i) = ra(k+2,i) - ddqk(i,k)*(epsln - ra(k,i)) ra(k,i) = epsln endif enddo enddo c !----------------------------------------------------------------------- ! at the kmax-1 level negative mixing ratios can only be borrowed ! from the kmax level !----------------------------------------------------------------------- do i = i1,i2 if(ra(kmax-1,i) .LE. epsln) then ra(kmax,i) = ra(kmax,i) - ddqk(i,kmax-1)* & (epsln - ra(kmax-1,i)) ra(kmax-1,i) = epsln endif if(ra(kmax,i) .LE. epsln) then ra(kmax,i) = epsln endif enddo if(toff .EQ.1.AND. (I1.EQ.1.OR.I2.EQ.imax))return if(toff .EQ.1.AND. (j.LE. 2+ngr.OR.J.GE.jmax-1-ngr))return if(itr .EQ. 0 .AND. l .EQ. 0) return !----------------------------------------------------------------------- ! estimate tendency of pstar*t and pstar*r due to adjustment ! store local adusted t & r into forecast field !----------------------------------------------------------------------- c c set up for radiation c epso = .01e-3 do i = i1 , i2 do k = 1 , kmax convt = tdif(k,i) phys1 = rf(k,i,2) phys2 = rf(k,i,3) phys3 = rf(k,i,4) if(convt .GT. 0.0 .OR. phys1 .GT. epso * .OR. phys2 .GT. epso .OR. phys3 .GT. epso) then tdt3(k,i) = 100.0 else tdt3(k,i) = 0.0 endif enddo enddo c do i = i1,i2 do k = 1,kmax ccc tdt3(k,i) = (ta(k,i) - told(k,i))*psf(i)*rdeltx tdr3(k,i) = (ra(k,i) - rold(k,i))*psf(i)*rdeltx if(ngd .EQ. 1) then stor1(k,i,j,14) = 86400.*(ta(k,i) - tf(k,i))/deltat stor1(k,i,j,15) = 86400.* tdif(k,i)/deltat endiF if(ngd .EQ. 2) then stor2(k,i,j,14) = 86400.*(ta(k,i) - tf(k,i))/deltat stor2(k,i,j,15) = 86400.* tdif(k,i)/deltat endif if(ngd .EQ. 3) then stor3(k,i,j,14) = 86400.*(ta(k,i) - tf(k,i))/deltat stor3(k,i,j,15) = 86400.* tdif(k,i)/deltat endif tf(k,i) = ta(k,i) rf(k,i,1) = ra(k,i) enddo enddo c print *,'itr in bulkajustEND=', itr return end subroutine bulkadjust