C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/fizhi/fizhi_moist.F,v 1.5 2004/07/07 20:07:47 molod Exp $ C $Name: $ #include "CPP_OPTIONS.h" subroutine moistio (ndmoist,istrip,npcs, . lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup, . pz,pl,ple,dpres,pkht,pkl,tz,qz,bi,bj,ntracer,ptracer, . qqz,dumoist,dvmoist,dtmoist,dqmoist, . im,jm,lm,ptop, . iras,rainlsp,rainconv,snowfall, . nswcld,cldtot_sw,cldras_sw,cldlsp_sw,nswlz,swlz, . nlwcld,cldtot_lw,cldras_lw,cldlsp_lw,nlwlz,lwlz, . lpnt,myid) #ifdef ALLOW_DIAGNOSTICS #include "diagnostics.h" #endif c Input Variables c --------------- integer ndmoist,istrip,npcs integer lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup real pz(im,jm),pl(im,jm,lm),ple(im,jm,lm+1),dpres(im,jm,lm) real pkht(im,jm,lm+1),pkl(im,jm,lm) real tz(im,jm,lm),qz(im,jm,lm,ntracer) integer bi,bj,ntracer,ptracer real qqz(im,jm,lm) real dumoist(im,jm,lm),dvmoist(im,jm,lm) real dtmoist(im,jm,lm),dqmoist(im,jm,lm,ntracer) integer im,jm,lm real ptop integer iras real rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm) integer nswcld,nswlz real cldlsp_sw(im,jm,lm),cldras_sw(im,jm,lm) real cldtot_sw(im,jm,lm),swlz(im,jm,lm) integer nlwcld,nlwlz real cldlsp_lw(im,jm,lm),cldras_lw(im,jm,lm) real cldtot_lw(im,jm,lm),lwlz(im,jm,lm) logical lpnt integer myid c Local Variables c --------------- integer ncrnd,nsecf real fracqq, rh,temp1,temp2,dum integer snowcrit parameter (fracqq = 0.1) real cldsr(im,jm,lm) real srcld(istrip,lm) real plev real cldnow,cldlsp_mem,cldras_mem,cldras,watnow,watmin,cldmin real cldprs(im,jm),cldtmp(im,jm) real cldhi (im,jm),cldlow(im,jm) real cldmid(im,jm),totcld(im,jm) real CLDLS(im,jm,lm) , CPEN(im,jm,lm) real tmpimjm(im,jm) real lsp_new(im,jm) real conv_new(im,jm) real snow_new(im,jm) real qqcolmin(im,jm) real qqcolmax(im,jm) integer levpbl(im,jm) c Gathered Arrays for Variable Cloud Base c --------------------------------------- real raincgath(im*jm) real pigather(im*jm) real thgather(im*jm,lm) real shgather(im*jm,lm) real pkzgather(im*jm,lm) real pkegather(im*jm,lm) real tmpgather(im*jm,lm) real deltgather(im*jm,lm) real delqgather(im*jm,lm) real ugather(im*jm,lm,ntracer) real delugather(im*jm,lm,ntracer) real deltrnev(im*jm,lm) real delqrnev(im*jm,lm) integer nindeces(lm) integer pblindex(im*jm) integer levgather(im*jm) c Stripped Arrays c --------------- real saveth (istrip,lm) real saveq (istrip,lm) real saveu (istrip,lm,ntracer) real usubcl (istrip, ntracer) real ple(istrip,lm+1), gam(istrip,lm) real TL(ISTRIP,lm) , SHL(ISTRIP,lm) real PL(ISTRIP,lm) , PLK(ISTRIP,lm) real PLKE(ISTRIP,lm+1) real TH(ISTRIP,lm) ,CVTH(ISTRIP,lm) real SHSAT(ISTRIP,lm) , CVQ(ISTRIP,lm) real UL(ISTRIP,lm,ntracer) real cvu(istrip,lm,ntracer) real CLMAXO(ISTRIP,lm),CLBOTH(ISTRIP,lm) real CLSBTH(ISTRIP,lm) real TMP1(ISTRIP,lm), TMP2(ISTRIP,lm) real TMP3(ISTRIP,lm), TMP4(ISTRIP,lm+1) real TMP5(ISTRIP,lm+1) integer ITMP1(ISTRIP,lm), ITMP2(ISTRIP,lm) integer ITMP3(ISTRIP,lm) real PRECIP(ISTRIP), PCMID(ISTRIP), PCNET(ISTRIP) real PCLOW (ISTRIP), SP(ISTRIP), PREP(ISTRIP) real PCPEN (ISTRIP,lm) integer pbl(istrip),depths(lm) real cldlz(istrip,lm), cldwater(im,jm,lm) real rhfrac(istrip), rhmin, pup, ppbl, rhcrit(istrip,lm) real offset, alpha, rasmax logical first logical lras real clfrac (istrip,lm) real cldmas (istrip,lm) real detrain(istrip,lm) real psubcld (istrip), psubcldg (im,jm) real psubcld_cnt(istrip), psubcldgc(im,jm) real rnd(lm/2) DATA FIRST /.TRUE./ integer imstp,nsubcl,nlras integer i,j,iloop,index,l,nn,num,numdeps,nt real tmstp,tminv,sday,grav,alhl,cp,elocp,gamfac real rkappa,p0kappa,p0kinv,ptopkap,pcheck real tice,getcon,pi C ********************************************************************** C **** INITIALIZATION **** C ********************************************************************** IMSTP = nsecf(NDMOIST) TMSTP = FLOAT(IMSTP) TMINV = 1. / TMSTP C Minimum Large-Scale Cloud Fraction at rhcrit alpha = 0.80 C Difference in fraction between SR and LS Threshold offset = 0.10 C Large-Scale Relative Humidity Threshold in PBL rhmin = 0.90 C Maximum Cloud Fraction associated with RAS rasmax = 1.00 nn = 3*3600.0/tmstp + 1 C Threshold for Cloud Fraction Memory cldmin = rasmax*(1.0-tmstp/3600.)**nn C Threshold for Cloud Liquid Water Memory watmin = 1.0e-8 SDAY = GETCON('SDAY') GRAV = GETCON('GRAVITY') ALHL = GETCON('LATENT HEAT COND') CP = GETCON('CP') ELOCP = GETCON('LATENT HEAT COND') / GETCON('CP') GAMFAC = GETCON('LATENT HEAT COND') * GETCON('EPS') * ELOCP . / GETCON('RGAS') RKAPPA = GETCON('KAPPA') P0KAPPA = 1000.0**RKAPPA P0KINV = 1. / P0KAPPA PTOPKAP = PTOP**RKAPPA tice = getcon('FREEZING-POINT') PI = 4.*atan(1.) c Determine Total number of Random Clouds to Check c --------------------------------------------- ncrnd = (lm-nltop+1)/2 if(first .and. myid.eq.0) then print * print *,'Top Level Allowed for Convection : ',nltop print *,' Highest Sub-Cloud Level: ',nsubmax print *,' Lowest Sub-Cloud Level: ',nsubmin print *,' Total Number of Random Clouds: ',ncrnd print * first = .false. endif c And now find PBL depth - the level where qq = fracqq * qq at surface c -------------------------------------------------------------------- do j = 1,jm do i = 1,im qqcolmin(i,j) = qqz(i,j,lm)*fracqq qqcolmax(i,j) = qqz(i,j,lm) levpbl(i,j) = lm+1 enddo enddo DO L = lm-1,1,-1 DO j = 1,jm DO i = 1,im IF((qqz(i,j,l).gt.qqcolmax(i,j)) 1 .and.(levpbl(i,j).eq.lm+1))then qqcolmax(i,j) = qqz(i,j,l) qqcolmin(i,j) = fracqq*qqcolmax(i,j) endif if((qqz(i,j,l).lt.qqcolmin(i,j)) 1 .and.(levpbl(i,j).eq.lm+1))levpbl(i,j)=L+1 enddo enddo enddo do j = 1,jm do i = 1,im if(levpbl(i,j).gt.nsubmin) levpbl(i,j) = nsubmin if(levpbl(i,j).lt.nsubmax) levpbl(i,j) = nsubmax enddo enddo c Set up the array of indeces of subcloud levels for the gathering c ---------------------------------------------------------------- index = 0 do L = nsubmin,nltop,-1 do j = 1,jm do i = 1,im if(levpbl(i,j).eq.L) then index = index + 1 pblindex(index) = (j-1)*im + i endif enddo enddo enddo do index = 1,im*jm levgather(index) = levpbl(pblindex(index),1) pigather(index) = pz(pblindex(index),1) enddo do L = 1,lm do index = 1,im*jm thgather(index,L) = tz(pblindex(index),1,L) shgather(index,L) = qz(pblindex(index),1,L,1) pkegather(index,L) = pkht(pblindex(index),1,L) pkzgather(index,L) = pkl (pblindex(index),1,L) enddo enddo do nt = 1,ntracer-ptracer do L = 1,lm do index = 1,im*jm ugather(index,L,nt) = qz(pblindex(index),1,L,nt+ptracer) enddo enddo enddo c bump the counter for number of calls to convection c -------------------------------------------------- iras = iras + 1 if( iras.ge.1e9 ) iras = 1 c select the 'random' cloud detrainment levels for RAS c ---------------------------------------------------- call rndcloud(iras,ncrnd,rnd,myid) do l=1,lm do j=1,jm do i=1,im dtmoist(i,j,l) = 0. do nt = 1,ntracer dqmoist(i,j,l,nt) = 0. enddo enddo enddo enddo C*********************************************************************** C **** LOOP OVER NPCS PEICES **** C ********************************************************************** DO 1000 NN = 1,NPCS C ********************************************************************** C **** VARIABLE INITIALIZATION **** C ********************************************************************** CALL STRIP ( pigather, SP ,im*jm,ISTRIP,1 ,NN ) CALL STRIP ( pkzgather, PLK ,im*jm,ISTRIP,lm,NN ) CALL STRIP ( pkegather, PLKE ,im*jm,ISTRIP,lm,NN ) CALL STRIP ( thgather, TH ,im*jm,ISTRIP,lm,NN ) CALL STRIP ( shgather, SHL ,im*jm,ISTRIP,lm,NN ) CALL STRINT( levgather, pbl ,im*jm,ISTRIP,1 ,NN ) do nt = 1,ntracer-ptracer call strip ( ugather(1,1,nt), ul(1,1,nt),im*jm,istrip,lm,nn ) enddo do l = 1,lm do i = 1,istrip PL(I,L) = SIG(L)*SP(I) + PTOP PLE(I,L) = SIGE(L)*SP(I) + PTOP enddo enddo do i = 1,istrip PLE(I,lm+1) = SP(I) + PTOP enddo C ********************************************************************** C **** SETUP FOR RAS CUMULUS PARAMETERIZATION **** C ********************************************************************** DO L = 1,lm DO I = 1,ISTRIP TH(I,L) = TH(I,L) * P0KAPPA CLMAXO(I,L) = 0. CLBOTH(I,L) = 0. cldmas(I,L) = 0. detrain(I,L) = 0. ENDDO ENDDO do L = 1,lm depths(L) = 0 enddo numdeps = 0 do L = nsubmin,nltop,-1 nindeces(L) = 0 do i = 1,istrip if(pbl(i).eq.L) nindeces(L) = nindeces(L) + 1 enddo if(nindeces(L).gt.0) then numdeps = numdeps + 1 depths(numdeps) = L endif enddo C Initiate a do-loop around RAS for the number of different C sub-cloud layer depths in this strip C --If all subcloud depths are the same, execute loop once C Otherwise loop over different subcloud layer depths num = 1 DO iloop = 1,numdeps nsubcl = depths(iloop) c Compute sub-cloud values for Temperature and Spec.Hum. c ------------------------------------------------------ DO 600 I=num,num+nindeces(nsubcl)-1 TMP1(I,2) = 0. TMP1(I,3) = 0. 600 CONTINUE NLRAS = NSUBCL - NLTOP + 1 DO 601 L=NSUBCL,lm DO 602 I=num,num+nindeces(nsubcl)-1 TMP1(I,2) = TMP1(I,2) + (PLE(I,L+1)-PLE(I,L))*TH (I,L)/sp(i) TMP1(I,3) = TMP1(I,3) + (PLE(I,L+1)-PLE(I,L))*SHL(I,L)/sp(i) 602 CONTINUE 601 CONTINUE DO 603 I=num,num+nindeces(nsubcl)-1 TMP1(I,4) = 1. / ( (PLE(I,lm+1)-PLE(I,NSUBCL))/sp(I) ) TH(I,NSUBCL) = TMP1(I,2)*TMP1(I,4) SHL(I,NSUBCL) = TMP1(I,3)*TMP1(I,4) 603 CONTINUE c Save initial value of tracers and compute sub-cloud value c --------------------------------------------------------- DO NT = 1,ntracer-ptracer do L = 1,lm do i = num,num+nindeces(nsubcl)-1 saveu(i,L,nt) = ul(i,L,nt) enddo enddo DO I=num,num+nindeces(nsubcl)-1 TMP1(I,2) = 0. ENDDO DO L=NSUBCL,lm DO I=num,num+nindeces(nsubcl)-1 TMP1(I,2) = TMP1(I,2)+(PLE(I,L+1)-PLE(I,L))*UL(I,L,NT)/sp(i) ENDDO ENDDO DO I=num,num+nindeces(nsubcl)-1 UL(I,NSUBCL,NT) = TMP1(I,2)*TMP1(I,4) usubcl(i,nt) = ul(i,nsubcl,nt) ENDDO ENDDO c Compute Pressure Arrays for RAS c ------------------------------- DO 111 L=1,lm DO 112 I=num,num+nindeces(nsubcl)-1 TMP4(I,L) = PLE(I,L) 112 CONTINUE 111 CONTINUE DO I=num,num+nindeces(nsubcl)-1 TMP5(I,1) = PTOPKAP / P0KAPPA ENDDO DO L=2,lm DO I=num,num+nindeces(nsubcl)-1 TMP5(I,L) = PLKE(I,L-1)*P0KINV ENDDO ENDDO DO I=num,num+nindeces(nsubcl)-1 TMP4(I,lm+1) = PLE (I,lm+1) TMP5(I,lm+1) = PLKE(I,lm)*P0KINV ENDDO DO 113 I=num,num+nindeces(nsubcl)-1 TMP4(I,NSUBCL+1) = PLE (I,lm+1) TMP5(I,NSUBCL+1) = PLKE(I,lm)*P0KINV 113 CONTINUE do i=num,num+nindeces(nsubcl)-1 C Temperature at top of sub-cloud layer tmp2(i,1) = TH(i,NSUBCL) * PLKE(i,NSUBCL)/P0KAPPA C Pressure at top of sub-cloud layer tmp2(i,2) = tmp4(i,nsubcl) enddo C CHANGED THIS: no RH requirement for RAS c call vqsat ( tmp2(num,1),tmp2(num,2),tmp2(num,3), c . dum,.false.,nindeces(nsubcl) ) c do i=num,num+nindeces(nsubcl)-1 c rh = SHL(I,NSUBCL) / tmp2(i,3) c if (rh .le. 0.85) then c rhfrac(i) = 0. c else if (rh .ge. 0.95) then c rhfrac(i) = 1. c else c rhfrac(i) = (rh-0.85)*10. c endif c enddo do i=num,num+nindeces(nsubcl)-1 rhfrac(i) = 1. enddo C Compute RH threshold for Large-scale condensation C Used in Slingo-Ritter clouds as well - define offset between SR and LS C Top level of atan func above this rh_threshold = rhmin pup = 600. do i=num,num+nindeces(nsubcl)-1 do L = nsubcl, lm rhcrit(i,L) = 1. enddo do L = 1, nsubcl-1 pcheck = (1000.-ptop)*sig(L) + ptop if (pcheck .le. pup) then rhcrit(i,L) = rhmin else ppbl = (1000.-ptop)*sig(nsubcl) + ptop rhcrit(i,L) = rhmin + (1.-rhmin)/(19.) * . ((atan( (2.*(pcheck-pup)/(ppbl-pup)-1.) * . tan(20.*pi/21.-0.5*pi) ) . + 0.5*pi) * 21./pi - 1.) endif enddo enddo c Save Initial Values of Temperature and Specific Humidity c -------------------------------------------------------- do L = 1,lm do i = num,num+nindeces(nsubcl)-1 saveth(i,L) = th (i,L) saveq (i,L) = shl(i,L) PCPEN (i,L) = 0. CLFRAC(i,L) = 0. enddo enddo CALL RAS ( NN,istrip,nindeces(nsubcl),NLRAS,NLTOP,lm,TMSTP 1, UL(num,1,1),ntracer-ptracer,TH(num,NLTOP),SHL(num,NLTOP) 2, TMP4(num,NLTOP), TMP5(num,NLTOP),rnd, ncrnd, PCPEN(num,NLTOP) 3, CLBOTH(num,NLTOP), CLFRAC(num,NLTOP) 4, cldmas(num,nltop), detrain(num,nltop) 8, cp,grav,rkappa,alhl,rhfrac(num),rasmax ) c Compute Diagnostic CLDMAS in RAS Subcloud Layers c ------------------------------------------------ do L=nsubcl,lm dum = dsig(L)/(1.0-sige(nsubcl)) do I=num,num+nindeces(nsubcl)-1 cldmas(i,L) = cldmas(i,L-1) - dum*cldmas(i,nsubcl-1) enddo enddo c Update Theta and Moisture due to RAS c ------------------------------------ DO L=1,nsubcl DO I=num,num+nindeces(nsubcl)-1 CVTH(I,L) = (TH (I,L) - saveth(i,l)) CVQ (I,L) = (SHL(I,L) - saveq (i,l)) ENDDO ENDDO DO L=nsubcl+1,lm DO I=num,num+nindeces(nsubcl)-1 CVTH(I,L) = cvth(i,nsubcl) CVQ (I,L) = cvq (i,nsubcl) ENDDO ENDDO DO L=nsubcl+1,lm DO I=num,num+nindeces(nsubcl)-1 TH (I,L) = saveth(i,l) + cvth(i,l) SHL(I,L) = saveq (i,l) + cvq (i,l) ENDDO ENDDO DO L=1,lm DO I=num,num+nindeces(nsubcl)-1 CVTH(I,L) = CVTH(I,L) *P0KINV*SP(I)*tminv CVQ (I,L) = CVQ (I,L) *SP(I)*tminv ENDDO ENDDO c Compute Tracer Tendency due to RAS c ---------------------------------- do nt = 1,ntracer-ptracer DO L=1,nsubcl-1 DO I=num,num+nindeces(nsubcl)-1 CVU(I,L,nt) = ( UL(I,L,nt)-saveu(i,l,nt) )*sp(i)*tminv ENDDO ENDDO DO L=nsubcl,lm DO I=num,num+nindeces(nsubcl)-1 if( usubcl(i,nt).ne.0.0 ) then cvu(i,L,nt) = ( ul(i,nsubcl,nt)-usubcl(i,nt) ) * . ( saveu(i,L,nt)/usubcl(i,nt) )*sp(i)*tminv else cvu(i,L,nt) = 0.0 endif ENDDO ENDDO enddo c Compute Diagnostic PSUBCLD (Subcloud Layer Pressure) c ---------------------------------------------------- do i=num,num+nindeces(nsubcl)-1 lras = .false. do L=nltop,nsubcl if( cvq(i,L).ne.0.0 ) lras = .true. enddo psubcld (i) = 0.0 psubcld_cnt(i) = 0.0 if( lras ) then psubcld (i) = sp(i)+ptop-ple(i,nsubcl) psubcld_cnt(i) = 1.0 endif enddo C End of subcloud layer depth loop (iloop) num = num+nindeces(nsubcl) ENDDO C ********************************************************************** C **** TENDENCY UPDATES **** C **** (Keep 'Gathered' tendencies in 'gather' arrays now) **** C ********************************************************************** call paste( CVTH,deltgather,istrip,im*jm,lm,NN ) call paste( CVQ,delqgather,istrip,im*jm,lm,NN ) do nt = 1,ntracer-ptracer call paste( CVU(1,1,nt),delugather(1,1,nt),istrip,im*jm,lm,NN ) enddo C ********************************************************************** C And now paste some arrays for filling diagnostics C (use pkegather to hold detrainment and tmpgather for cloud mass flux) C ********************************************************************** if(icldmas .gt.0) call paste( cldmas,tmpgather,istrip,im*jm,lm,NN) if(idtrain .gt.0) call paste(detrain,pkegather,istrip,im*jm,lm,NN) if(ipsubcld.gt.0) then call paste(psubcld ,psubcldg ,istrip,im*jm,1,NN) call paste(psubcld_cnt,psubcldgc,istrip,im*jm,1,NN) endif C ********************************************************************* C **** RE-EVAPORATION OF PENETRATING CONVECTIVE RAIN **** C ********************************************************************* CALL STRIP ( thgather,TH ,im*jm,ISTRIP,lm,NN) CALL STRIP ( shgather,SHL,im*jm,ISTRIP,lm,NN) DO L=1,lm DO I=1,ISTRIP TH(I,L) = TH(I,L) + CVTH(I,L)*tmstp/SP(I) SHL(I,L) = SHL(I,L) + CVQ(I,L)*tmstp/SP(I) TL(I,L) = TH(I,L)*PLK(I,L) saveth(I,L) = th(I,L) saveq (I,L) = SHL(I,L) ENDDO ENDDO CALL RNEVP (NN,ISTRIP,lm,TL,SHL,PCPEN,PL,CLFRAC,SP,DSIG,PLKE, . PLK,TH,TMP1,TMP2,TMP3,ITMP1,ITMP2,PCNET,PRECIP, . CLSBTH,TMSTP,1.,cp,grav,alhl,gamfac,cldlz,rhcrit,offset,alpha) C ********************************************************************** C **** TENDENCY UPDATES **** C ********************************************************************** DO L=1,lm DO I =1,ISTRIP TMP1(I,L) = sp(i) * (SHL(I,L)-saveq(I,L)) * tminv ENDDO CALL PSTBMP(TMP1(1,L),delqgather(1,L),ISTRIP,im*jm,1,NN) DO I =1,ISTRIP TMP1(I,L) = sp(i) * ((TL(I,L)/PLK(I,L))-saveth(i,l)) * tminv ENDDO CALL PSTBMP(TMP1(1,L),deltgather(1,L),ISTRIP,im*jm,1,NN) C Paste rain evap tendencies into arrays for diagnostic output c ------------------------------------------------------------ if(idtls.gt.0)then DO I =1,ISTRIP TMP1(I,L) = ((TL(I,L)/PLK(I,L))-saveth(i,l))*plk(i,l)*sday*tminv ENDDO call paste(tmp1(1,L),deltrnev(1,L),istrip,im*jm,1,NN) endif if(idqls.gt.0)then DO I =1,ISTRIP TMP1(I,L) = (SHL(I,L)-saveq(I,L)) * 1000. * sday * tminv ENDDO call paste(tmp1(1,L),delqrnev(1,L),istrip,im*jm,1,NN) endif ENDDO C ********************************************************************* C Add Non-Precipitating Clouds where the relative C humidity is less than 100% C Apply Cloud Top Entrainment Instability C ********************************************************************* do L=1,lm do i=1,istrip srcld(i,L) = -clsbth(i,L) enddo enddo call srclouds (saveth,saveq,plk,pl,plke,clsbth,cldlz,istrip,lm, . rhcrit,offset,alpha) do L=1,lm do i=1,istrip srcld(i,L) = srcld(i,L)+clsbth(i,L) enddo enddo C ********************************************************************* C **** PASTE CLOUD AMOUNTS **** C ********************************************************************* call paste ( srcld, cldsr,istrip,im*jm,lm,nn ) call paste ( cldlz,cldwater,istrip,im*jm,lm,nn ) call paste ( clsbth, cldls,istrip,im*jm,lm,nn ) call paste ( clboth, cpen ,istrip,im*jm,lm,nn ) c compute Total Accumulated Precip for Landsurface Model c ------------------------------------------------------ do i = 1,istrip C Initialize Rainlsp, Rainconv and Snowfall tmp1(i,1) = 0.0 tmp1(i,2) = 0.0 tmp1(i,3) = 0.0 enddo do i = 1,istrip prep(i) = PRECIP(I) + PCNET(I) tmp1(i,1) = PRECIP(I) tmp1(i,2) = pcnet(i) enddo c c check whether there is snow c------------------------------------------------------- c snow algorthm: c if temperature profile from the surface level to 700 mb c uniformaly c below zero, then precipitation (total) is c snowfall. Else there is no snow. c------------------------------------------------------- do i = 1,istrip snowcrit=0 do l=lup,lm if (saveth(i,l)*plk(i,l).le. tice ) then snowcrit=snowcrit+1 endif enddo if (snowcrit .eq. (lm-lup+1)) then tmp1(i,3) = prep(i) tmp1(i,1)=0.0 tmp1(i,2)=0.0 endif enddo CALL paste (tmp1(1,1), lsp_new,ISTRIP,im*jm,1,NN) CALL paste (tmp1(1,2),conv_new,ISTRIP,im*jm,1,NN) CALL paste (tmp1(1,3),snow_new,ISTRIP,im*jm,1,NN) if(iprecon.gt.0) then CALL paste (pcnet,raincgath,ISTRIP,im*jm,1,NN) endif C ********************************************************************* C **** End Major Stripped Region **** C ********************************************************************* 1000 CONTINUE C Large Scale Rainfall, Conv rain, and snowfall c --------------------------------------------- call back2grd ( lsp_new,pblindex, lsp_new,im*jm) call back2grd (conv_new,pblindex,conv_new,im*jm) call back2grd (snow_new,pblindex,snow_new,im*jm) if(iprecon.gt.0) then call back2grd (raincgath,pblindex,raincgath,im*jm) endif c Subcloud Layer Pressure c ----------------------- if(ipsubcld.gt.0) then call back2grd (psubcldg ,pblindex,psubcldg ,im*jm) call back2grd (psubcldgc,pblindex,psubcldgc,im*jm) endif do L = 1,lm C Delta theta,q, convective, max and ls clouds c -------------------------------------------- call back2grd (deltgather(1,L),pblindex, dtmoist(1,1,L) ,im*jm) call back2grd (delqgather(1,L),pblindex, dqmoist(1,1,L,1),im*jm) call back2grd ( cpen(1,1,L),pblindex, cpen(1,1,L) ,im*jm) call back2grd ( cldls(1,1,L),pblindex, cldls(1,1,L) ,im*jm) call back2grd (cldwater(1,1,L),pblindex,cldwater(1,1,L) ,im*jm) call back2grd ( pkzgather(1,L),pblindex, pkzgather(1,L) ,im*jm) C Diagnostics: c ------------ if(icldmas.gt.0)call back2grd(tmpgather(1,L),pblindex, . tmpgather(1,L),im*jm) if(idtrain.gt.0)call back2grd(pkegather(1,L),pblindex, . pkegather(1,L),im*jm) if(idtls.gt.0)call back2grd(deltrnev(1,L),pblindex, . deltrnev(1,L),im*jm) if(idqls.gt.0)call back2grd(delqrnev(1,L),pblindex, . delqrnev(1,L),im*jm) if(icldnp.gt.0)call back2grd(cldsr(1,1,L),pblindex, . cldsr(1,1,L),im*jm) enddo c Tracers c ------- do nt = 1,ntracer-ptracer do L = 1,lm call back2grd (delugather(1,L,nt),pblindex, . dqmoist(1,1,L,ptracer+nt),im*jm) enddo enddo C ********************************************************************** C BUMP DIAGNOSTICS C ********************************************************************** c Clear-Sky (Above 400mb) Temperature c ----------------------------------- if( itmpuclr.ne.0 .or. isphuclr.ne.0 ) then do j = 1,jm do i = 1,im totcld(i,j) = 0.0 enddo enddo do L = 1,midlevel do j = 1,jm do i = 1,im if(cldls(i,j,L).ne.0.0.or.cpen(i,j,L).ne.0.0)totcld(i,j) = 1.0 enddo enddo enddo do L = 1,lm if( itmpuclr.ne.0 ) then do i = 1,im*jm if( totcld(i,1).eq.0.0 ) then qdiag(i,1,itmpuclr +L-1,bi,bj) = . qdiag(i,1,itmpuclr +L-1,bi,bj) + tz(i,1,L)*pkzgather(i,L) qdiag(i,1,itmpuclrc+L-1,bi,bj) = . qdiag(i,1,itmpuclrc+L-1,bi,bj)+1.0 endif enddo endif if( isphuclr.ne.0 ) then do i = 1,im*jm if( totcld(i,1).eq.0.0 ) then qdiag(i,1,isphuclr +L-1,bi,bj) = . qdiag(i,1,isphuclr +L-1,bi,bj) + qz(i,1,L,1)*1000.0 qdiag(i,1,isphuclrc+L-1,bi,bj) = . qdiag(i,1,isphuclrc+L-1,bi,bj) + 1.0 endif enddo endif enddo endif c Sub-Cloud Layer c ------------------------- if( ipsubcld.ne.0 ) then do j = 1,jm do i = 1,im qdiag(i,j,ipsubcld,bi,bj) = qdiag(i,j,ipsubcld,bi,bj) + . psubcldg (i,j) qdiag(i,j,ipsubcldc,bi,bj) = qdiag(i,j,ipsubcldc,bi,bj) + . psubcldgc(i,j) enddo enddo endif c Non-Precipitating Cloud Fraction c -------------------------------- if( icldnp.ne.0 ) then do L = 1,lm do j = 1,jm do i = 1,im qdiag(i,j,icldnp+L-1,bi,bj) = qdiag(i,j,icldnp+L-1,bi,bj) + . cldsr(i,j,L) enddo enddo enddo ncldnp = ncldnp + 1 endif c Moist Processes Heating Rate c ---------------------------- if(imoistt.gt.0) then do L = 1,lm do i = 1,im*jm qdiag(i,1,imoistt+L-1,bi,bj) = qdiag(i,1,imoistt+L-1,bi,bj) + . (dtmoist(i,1,L)*sday*pkzgather(i,L)/pz(i,1)) enddo enddo endif c Moist Processes Moistening Rate c ------------------------------- if(imoistq.gt.0) then do L = 1,lm do j = 1,jm do i = 1,im qdiag(i,j,imoistq+L-1,bi,bj) = qdiag(i,j,imoistq+L-1,bi,bj) + . (dqmoist(i,j,L,1)*sday*1000.0/pz(i,j)) enddo enddo enddo endif c Cloud Mass Flux c --------------- if(icldmas.gt.0) then do L = 1,lm do i = 1,im*jm qdiag(i,1,icldmas+L-1,bi,bj) = qdiag(i,1,icldmas+L-1,bi,bj) + . tmpgather(i,L) enddo enddo endif c Detrained Cloud Mass Flux c ------------------------- if(idtrain.gt.0) then do L = 1,lm do i = 1,im*jm qdiag(i,1,idtrain+L-1,bi,bj) = qdiag(i,1,idtrain+L-1,bi,bj) + . pkegather(i,L) enddo enddo endif c Grid-Scale Condensational Heating Rate c -------------------------------------- if(idtls.gt.0) then do L = 1,lm do i = 1,im*jm qdiag(i,1,idtls+L-1,bi,bj) = qdiag(i,1,idtls+L-1,bi,bj) + . deltrnev(i,L) enddo enddo endif c Grid-Scale Condensational Moistening Rate c ----------------------------------------- if(idqls.gt.0) then do L = 1,lm do i = 1,im*jm qdiag(i,1,idqls+L-1,bi,bj) = qdiag(i,1,idqls+L-1,bi,bj) + . delqrnev(i,L) enddo enddo endif c Total Precipitation c ------------------- if(ipreacc.gt.0) then do j = 1,jm do i = 1,im qdiag(i,j,ipreacc,bi,bj) = qdiag(i,j,ipreacc,bi,bj) . + ( lsp_new(I,j) . + snow_new(I,j) . + conv_new(i,j) ) *sday*tminv enddo enddo endif c Convective Precipitation c ------------------------ if(iprecon.gt.0) then do i = 1,im*jm qdiag(i,1,iprecon,bi,bj) = qdiag(i,1,iprecon,bi,bj) + . raincgath(i)*sday*tminv enddo endif C ********************************************************************** C **** Fill Rainfall and Snowfall Arrays for Land Surface Model **** C **** Note: Precip Rates work when DT(turb) Cloud Fraction from RAS *** C *** CLDLS => Cloud Fraction from RNEVP *** C ********************************************************************** do j = 1,jm do i = 1,im cldhi (i,j) = 0. cldmid(i,j) = 0. cldlow(i,j) = 0. cldtmp(i,j) = 0. cldprs(i,j) = 0. tmpimjm(i,j) = 0. enddo enddo c Set Moist-Process Memory Coefficient c ------------------------------------ cldras_mem = 1.0-tmstp/ 3600.0 cldlsp_mem = 1.0-tmstp/(3600.0*3) do L = 1,lm do i = 1,im*jm plev = sig(L)*pz(i,1)+ptop c Compute Time-averaged Cloud and Water Amounts for Longwave Radiation c -------------------------------------------------------------------- watnow = cldwater(i,1,L) if( plev.le.500.0 ) then cldras = min( max( cldras_lw(i,1,L)*cldras_mem,cpen(i,1,L)),1.0) else cldras = 0.0 endif cldlsp = min( max( cldlsp_lw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0) if( cldras.lt.cldmin ) cldras = 0.0 if( cldlsp.lt.cldmin ) cldlsp = 0.0 cldnow = max( cldlsp,cldras ) lwlz(i,1,L) = ( nlwlz*lwlz(i,1,L) + watnow)/(nlwlz +1) cldtot_lw(i,1,L) = (nlwcld*cldtot_lw(i,1,L) + cldnow)/(nlwcld+1) cldlsp_lw(i,1,L) = (nlwcld*cldlsp_lw(i,1,L) + cldlsp)/(nlwcld+1) cldras_lw(i,1,L) = (nlwcld*cldras_lw(i,1,L) + cldras)/(nlwcld+1) c Compute Time-averaged Cloud and Water Amounts for Shortwave Radiation c --------------------------------------------------------------------- watnow = cldwater(i,1,L) if( plev.le.500.0 ) then cldras = min( max(cldras_sw(i,1,L)*cldras_mem, cpen(i,1,L)),1.0) else cldras = 0.0 endif cldlsp = min( max(cldlsp_sw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0) if( cldras.lt.cldmin ) cldras = 0.0 if( cldlsp.lt.cldmin ) cldlsp = 0.0 cldnow = max( cldlsp,cldras ) swlz(i,1,L) = ( nswlz*swlz(i,1,L) + watnow)/(nswlz +1) cldtot_sw(i,1,L) = (nswcld*cldtot_sw(i,1,L) + cldnow)/(nswcld+1) cldlsp_sw(i,1,L) = (nswcld*cldlsp_sw(i,1,L) + cldlsp)/(nswcld+1) cldras_sw(i,1,L) = (nswcld*cldras_sw(i,1,L) + cldras)/(nswcld+1) c Compute Instantaneous Low-Mid-High Maximum Overlap Cloud Fractions c ---------------------------------------------------------------------- if( L.lt.midlevel ) cldhi (i,1) = max( cldnow,cldhi (i,1) ) if( L.ge.midlevel .and. . L.lt.lowlevel ) cldmid(i,1) = max( cldnow,cldmid(i,1) ) if( L.ge.lowlevel ) cldlow(i,1) = max( cldnow,cldlow(i,1) ) c Compute Cloud-Top Temperature and Pressure c ------------------------------------------ cldtmp(i,1) = cldtmp(i,1) + cldnow*pkzgather(i,L) . * ( tz(i,1,L) + dtmoist(i,1,L)*tmstp/pz(i,1) ) cldprs(i,1) = cldprs(i,1) + cldnow*plev tmpimjm(i,1) = tmpimjm(i,1) + cldnow enddo enddo c Compute Instantanious Total 2-D Cloud Fraction c ---------------------------------------------- do j = 1,jm do i = 1,im totcld(i,j) = 1.0 - (1.-cldhi (i,j)) . * (1.-cldmid(i,j)) . * (1.-cldlow(i,j)) enddo enddo C ********************************************************************** C *** Fill Cloud Top Pressure and Temperature Diagnostic *** C ********************************************************************** if(icldtmp.gt.0) then do j = 1,jm do i = 1,im if( cldtmp(i,j).gt.0.0 ) then qdiag(i,j,icldtmp,bi,bj) = qdiag(i,j,icldtmp,bi,bj) + . cldtmp(i,j)*totcld(i,j)/tmpimjm(i,j) qdiag(i,j,icttcnt,bi,bj) = qdiag(i,j,icttcnt,bi,bj) + . totcld(i,j) endif enddo enddo endif if(icldprs.gt.0) then do j = 1,jm do i = 1,im if( cldprs(i,j).gt.0.0 ) then qdiag(i,j,icldprs,bi,bj) = qdiag(i,j,icldprs,bi,bj) + . cldprs(i,j)*totcld(i,j)/tmpimjm(i,j) qdiag(i,j,ictpcnt,bi,bj) = qdiag(i,j,ictpcnt,bi,bj) + . totcld(i,j) endif enddo enddo endif C ********************************************************************** C **** INCREMENT COUNTERS **** C ********************************************************************** nlwlz = nlwlz + 1 nswlz = nswlz + 1 nlwcld = nlwcld + 1 nswcld = nswcld + 1 nmoistt = nmoistt + 1 nmoistq = nmoistq + 1 npreacc = npreacc + 1 nprecon = nprecon + 1 ncldmas = ncldmas + 1 ndtrain = ndtrain + 1 ndtls = ndtls + 1 ndqls = ndqls + 1 RETURN END SUBROUTINE RAS( NN, LEN, LENC, K, NLTOP, nlayr, DT *, UOI, ntracer, POI, QOI, PRS, PRJ, rnd, ncrnd *, RAINS, CLN, CLF, cldmas, detrain *, cp,grav,rkappa,alhl,rhfrac,rasmax ) C C********************************************************************* C*********************** ARIES MODEL ******************************* C********************* SUBROUTINE RAS ***************************** C********************** 16 MARCH 1988 ****************************** C********************************************************************* C PARAMETER (KRMIN=01) PARAMETER (ICM=1000) PARAMETER (CMB2PA=100.0) PARAMETER (rknob = 10.) C integer ntracer integer nltop,nlayr DIMENSION UOI(len,nlayr,ntracer), POI(len,K) DIMENSION QOI(len,K), PRS(len,K+1), PRJ(len,K+1) dimension rnd(ncrnd) C DIMENSION RAINS(len,K), CLN(len,K), CLF(len,K) DIMENSION cldmas(len,K), detrain(len,K) DIMENSION TCU(len,K), QCU(len,K) real ucu(len,K,ntracer) DIMENSION ALF(len,K), BET(len,K), GAM(len,K) *, ETA(len,K), HOI(len,K) *, PRH(len,K), PRI(len,K) DIMENSION HST(len,K), QOL(len,K), GMH(len,K) DIMENSION TX1(len), TX2(len), TX3(len), TX4(len), TX5(len) *, TX6(len), TX7(len), TX8(len), TX9(len) *, TX11(len), TX12(len), TX13(len), TX14(len,ntracer) *, TX15(len), TX16(len) *, WFN(len), IA1(len), IA2(len), IA3(len) DIMENSION cloudn(len), pcu(len) real rhfrac(len),rasmax DIMENSION IC(ICM), IRND(icm) dimension cmass(len,K) LOGICAL SETRAS do L = 1,k do I = 1,LENC rains(i,l) = 0. enddo enddo p00 = 1000. crtmsf = 0. C The numerator here is the fraction of the subcloud layer mass flux C allowed to entrain into the cloud CCC FRAC = 1./dt FRAC = 0.5/dt KM1 = K - 1 KP1 = K + 1 C we want the ras adjustment time scale to be one hour (indep of dt) RASBLF = 1./3600. C KPRV = KM1 C Removed KRMAX parameter KCR = MIN(KM1,nlayr-2) KFX = KM1 - KCR NCMX = KFX + NCRND C IF (KFX .GT. 0) THEN DO NC=1,KFX IC(NC) = K - NC ENDDO ENDIF C IF (NCRND .GT. 0) THEN DO I=1,ncrnd IRND(I) = (RND(I)-0.0005)*(KCR-KRMIN+1) IRND(I) = IRND(I) + KRMIN ENDDO C DO NC=1,NCRND IC(KFX+NC) = IRND(NC) ENDDO ENDIF C DO 100 NC=1,NCMX C IF (NC .EQ. 1 ) THEN SETRAS = .TRUE. ELSE SETRAS = .FALSE. ENDIF IB = IC(NC) c Initialize Cloud Fraction Array c ------------------------------- do i = 1,lenc cloudn(i) = 0.0 enddo CALL CLOUD(nn,LEN, LENC, K, NLTOP, nlayr, IB, RASBLF,SETRAS,FRAC *, CP, ALHL, RKAPPA, GRAV, P00, CRTMSF *, POI, QOI, UOI, Ntracer, PRS, PRJ *, PCU, CLOUDN, TCU, QCU, UCU, CMASS *, ALF, BET, GAM, PRH, PRI, HOI, ETA *, HST, QOL, GMH *, TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9 *, WFN, TX11, TX12, TX13, TX14, TX15 *, IA1,IA2,IA3,rhfrac) C Compute fraction of grid box into which rain re-evap occurs (clf) c ----------------------------------------------------------------- do i = 1,lenc c mass in detrainment layer c ------------------------- tx1(i) = cmb2pa * (prs(i,ib+1) - prs(i,ib))/(grav*dt) c ratio of detraining cloud mass to mass in detrainment layer c ----------------------------------------------------------- tx1(i) = rhfrac(i)*rknob * cmass(i,ib) / tx1(i) if(cmass(i,K).gt.0.) clf(i,ib) = clf(i,ib) + tx1(i) if( clf(i,ib).gt.1.) clf(i,ib) = 1. enddo c Compute Total Cloud Mass Flux c ***************************** do L=ib,k do i=1,lenc cmass(i,L) = rhfrac(i)*cmass(i,L) * dt enddo enddo do L=ib,k do i=1,lenc cldmas(i,L) = cldmas(i,L) + cmass(i,L) enddo enddo do i=1,lenc detrain(i,ib) = detrain(i,ib) + cmass(i,ib) enddo DO L=IB,K DO I=1,LENC POI(I,L) = POI(I,L) + TCU(I,L) * DT * rhfrac(i) QOI(I,L) = QOI(I,L) + QCU(I,L) * DT * rhfrac(i) ENDDO ENDDO DO NT=1,Ntracer DO L=IB,K DO I=1,LENC UOI(I,L+nltop-1,NT)=UOI(I,L+nltop-1,NT)+UCU(I,L,NT)*DT*rhfrac(i) ENDDO ENDDO ENDDO DO I=1,LENC rains(I,ib) = rains(I,ib) + PCU(I)*dt * rhfrac(i) ENDDO 100 CONTINUE c Fill Convective Cloud Fractions based on 3-D Rain Amounts c --------------------------------------------------------- do L=k-1,1,-1 do i=1,lenc tx1(i) = 100*(prs(i,L+1)-prs(i,L))/grav cln(i,L) = min(1600*rains(i,L)/tx1(i),rasmax ) enddo enddo RETURN END subroutine rndcloud (iras,nrnd,rnd,myid) implicit none integer n,iras,nrnd,myid real random_numbx real rnd(nrnd) integer irm parameter (irm = 1000) real random(irm) integer i,mcheck,numrand,iseed,index logical first data first /.true./ integer iras0 data iras0 /0/ save random, iras0 if(nrnd.eq.0.)then do i = 1,nrnd rnd(i) = 0 enddo if(first .and. myid.eq.0) print *,' NO RANDOM CLOUDS IN RAS ' go to 100 endif mcheck = mod(iras-1,irm/nrnd) c First Time In From a Continuing RESTART (IRAS.GT.1) or Reading a New RESTART c ---------------------------------------------------------------------------- if( first.and.(iras.gt.1) .or. iras.ne.iras0+1 )then if( myid.eq.0 ) print *, 'Recreating Rand Numb Array in RNDCLOUD' if( myid.eq.0 ) print *, 'IRAS: ',iras,' IRAS0: ',iras0 numrand = mod(iras,irm/nrnd) * nrnd iseed = iras * nrnd - numrand call random_seedx(iseed) do i = 1,irm random(i) = random_numbx() enddo index = (iras-1)*nrnd c Multiple Time In But have Used Up all 1000 numbers (MCHECK.EQ.0) c ---------------------------------------------------------------- else if (mcheck.eq.0) then iseed = (iras-1)*nrnd call random_seedx(iseed) do i = 1,irm random(i) = random_numbx() enddo index = iseed c Multiple Time In But have NOT Used Up all 1000 numbers (MCHECK.NE.0) c -------------------------------------------------------------------- else index = (iras-1)*nrnd endif index = mod(index,irm) if( index+nrnd.gt.1000 ) index=1000-nrnd do n = 1,nrnd rnd(n) = random(index+n) enddo 100 continue first = .false. iras0 = iras return end real function random_numbx() implicit none #if CRAY real ranf random_numbx = ranf() #endif #if SGI real rand random_numbx = rand() #endif return end subroutine random_seedx (iseed) implicit none integer iseed #if CRAY call ranset (iseed) #endif #if SGI integer*4 seed seed = iseed call srand (seed) #endif return end SUBROUTINE CLOUD(nn,LEN, LENC, K, NLTOP, nlayr, IC, RASALF, *, SETRAS, FRAC *, CP, ALHL, RKAP, GRAV, P00, CRTMSF *, POI, QOI, UOI, Ntracer, PRS, PRJ *, PCU, CLN, TCU, QCU, UCU, CMASS *, ALF, BET, GAM, PRH, PRI, HOL, ETA *, HST, QOL, GMH *, TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, ALM *, WFN, AKM, QS1, CLF, UHT, WLQ *, IA, I1, I2,rhfrac) C C********************************************************************* C******************** Relaxed Arakawa-Schubert *********************** C********************* Plug Compatible Version ********************** C************************ SUBROUTINE CLOUD *************************** C************************* 23 JULY 1992 *************************** C********************************************************************* C********************************************************************* C********************************************************************* C************************** Developed By ***************************** C************************** ***************************** C************************ Shrinivas Moorthi ************************** C************************ and ************************** C************************ Max J. Suarez ***************************** C************************ ***************************** C******************** Laboratory for Atmospheres ********************* C****************** NASA/GSFC, Greenbelt, MD 20771 ******************* C********************************************************************* C********************************************************************* C C The calculations of Moorthi and Suarez (1992, MWR) are C contained in the CLOUD routine. C It is probably advisable, at least initially, to treat CLOUD C as a black box that computes the single cloud adjustments. RAS, C on the other hand, can be tailored to each GCMs configuration C (ie, number and placement of levels, nature of boundary layer, C time step and frequency with which RAS is called). C C C Input: C ------ C C LEN : The inner dimension of update and input arrays. C C LENC : The run: the number of soundings processes in a single call. C RAS works on the first LENC of the LEN soundings C passed. This allows working on pieces of the world C say for multitasking, without declaring temporary arrays C and copying the data to and from them. This is an f77 C version. An F90 version would have to allow more C flexibility in the argument declarations. Obviously C (LENC<=LEN). C C K : Number of vertical layers (increasing downwards). C Need not be the same as the number of layers in the C GCM, since it is the outer dimension. The bottom layer C (K) is the subcloud layer. C C IC : Detrainment level to check for presence of convection C C RASALF : Relaxation parameter (< 1.) for present cloud-type C C SETRAS : Logical parameter to control re-calculation of C saturation specific humidity and mid level P**kappa C C FRAC : Fraction of the PBL (layer K) mass allowed to be used C by a cloud-type in time DT C C CP : Specific heat at constant pressure C C ALHL : Latent Heat of condensation C C RKAP : R/Cp, where R is the gas constant C C GRAV : Acceleration due to gravity C C P00 : A reference pressure in hPa, useually 1000 hPa C C CRTMSF : Critical value of mass flux above which cloudiness at C the detrainment layer of that cloud-type is assumed. C Affects only cloudiness calculation. C C POI : 2D array of dimension (LEN,K) containing potential C temperature. Updated but not initialized by RAS. C C QOI : 2D array of dimension (LEN,K) containing specific C humidity. Updated but not initialized by RAS. C C UOI : 3D array of dimension (LEN,K,NTRACER) containing tracers C Updated but not initialized by RAS. C C PRS : 2D array of dimension (LEN,K+1) containing pressure C in hPa at the interfaces of K-layers from top of the C atmosphere to the bottom. Not modified. C C PRJ : 2D array of dimension (LEN,K+1) containing (PRS/P00) ** C RKAP. i.e. Exner function at layer edges. Not modified. C C rhfrac : 1D array of dimension (LEN) containing a rel.hum. scaling C fraction. Not modified. C C Output: C ------- C C PCU : 1D array of length LEN containing accumulated C precipitation in mm/sec. C C CLN : 2D array of dimension (LEN,K) containing cloudiness C Note: CLN is bumped but NOT initialized C C TCU : 2D array of dimension (LEN,K) containing accumulated C convective heating (K/sec). C C QCU : 2D array of dimension (LEN,K) containing accumulated C convective drying (kg/kg/sec). C C CMASS : 2D array of dimension (LEN,K) containing the C cloud mass flux (kg/sec). Filled from cloud top C to base. C C Temporaries: C C ALF, BET, GAM, ETA, PRH, PRI, HOI, HST, QOL, GMH are temporary C 2D real arrays of dimension of at least (LENC,K) where LENC is C the horizontal dimension over which convection is invoked. C C C TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9, AKM, QS1, CLF, UHT C VHT, WLQ WFN are temporary real arrays of length at least LENC C C IA, I1, and I2 are temporary integer arrays of length LENC C C C************************************************************************ C C PARAMETER (DAYLEN=86400.0, HALF=0.5, ONE=1.0, ZERO=0.0) PARAMETER (CMB2PA=100.0) PARAMETER (RHMAX=0.9999) C integer nltop,ntracer,nlayr DIMENSION POI(LEN,K), QOI(LEN,K), PRS(LEN,K+1) *, PRJ(LEN,K+1) *, TCU(LEN,K), QCU(LEN,K), CMASS(LEN,K), CLN(LEN) real uoi(len,nlayr,ntracer) DIMENSION ALF(LEN,K), BET(LEN,K), GAM(LEN,K) *, PRH(LEN,K), PRI(LEN,K) DIMENSION AKM(LENC), WFN(LENC) DIMENSION HOL(LENC,K), QOL(LENC,K), ETA(LENC,K), HST(LENC,K) *, GMH(LENC,K), ALM(LENC), WLQ(LENC), QS1(LENC) *, TX1(LENC), TX2(LENC), TX3(LENC), TX4(LENC) *, TX5(LENC), TX6(LENC), TX7(LENC), TX8(LENC) *, CLF(LENC), PCU(LENC) DIMENSION IA(LENC), I1(LENC), I2(LENC) real rhfrac(len) real ucu(len,k,ntracer),uht(len,ntracer) LOGICAL SETRAS integer nt c Explicit Inline Directives c -------------------------- #if CRAY #if f77 cfpp$ expand (qsat) #endif #endif RKAPP1 = 1.0 + RKAP ONEBCP = 1.0 / CP ALBCP = ALHL * ONEBCP ONEBG = 1.0 / GRAV CPBG = CP * ONEBG TWOBAL = 2.0 / ALHL C KM1 = K - 1 IC1 = IC + 1 C C SETTIING ALF, BET, GAM, PRH, AND PRI : DONE ONLY WHEN SETRAS=.T. C IF (SETRAS) THEN DO 2050 L=1,K DO 2030 I=1,LENC PRH(I,L) = (PRJ(I,L+1)*PRS(I,L+1) - PRJ(I,L)*PRS(I,L)) * / ((PRS(I,L+1)-PRS(I,L)) * RKAPP1) 2030 CONTINUE 2050 CONTINUE DO 2070 L=1,K DO 2060 I=1,LENC TX5(I) = POI(I,L) * PRH(I,L) TX1(I) = (PRS(I,L) + PRS(I,L+1)) * 0.5 TX3(I) = TX5(I) CALL QSAT(TX3(I), TX1(I), TX2(I), TX4(I), .TRUE.) ALF(I,L) = TX2(I) - TX4(I) * TX5(I) BET(I,L) = TX4(I) * PRH(I,L) GAM(I,L) = 1.0 / ((1.0 + TX4(I)*ALBCP) * PRH(I,L)) PRI(I,L) = (CP/CMB2PA) / (PRS(I,L+1) - PRS(I,L)) 2060 CONTINUE 2070 CONTINUE ENDIF C C DO 10 L=1,K DO 10 I=1,LEN TCU(I,L) = 0.0 QCU(I,L) = 0.0 CMASS(I,L) = 0.0 10 CONTINUE do nt = 1,ntracer do L=1,K do I=1,LENC ucu(I,L,nt) = 0.0 enddo enddo enddo C DO 30 I=1,LENC TX1(I) = PRJ(I,K+1) * POI(I,K) QS1(I) = ALF(I,K) + BET(I,K)*POI(I,K) QOL(I,K) = MIN(QS1(I)*RHMAX,QOI(I,K)) HOL(I,K) = TX1(I)*CP + QOL(I,K)*ALHL ETA(I,K) = ZERO TX2(I) = (PRJ(I,K+1) - PRJ(I,K)) * POI(I,K) * CP 30 CONTINUE C IF (IC .LT. KM1) THEN DO 3703 L=KM1,IC1,-1 DO 50 I=1,LENC QS1(I) = ALF(I,L) + BET(I,L)*POI(I,L) QOL(I,L) = MIN(QS1(I)*RHMAX,QOI(I,L)) C TEM1 = TX2(I) + PRJ(I,L+1) * POI(I,L) * CP HOL(I,L) = TEM1 + QOL(I,L )* ALHL HST(I,L) = TEM1 + QS1(I) * ALHL TX1(I) = (PRJ(I,L+1) - PRJ(I,L)) * POI(I,L) ETA(I,L) = ETA(I,L+1) + TX1(I)*CPBG TX2(I) = TX2(I) + TX1(I)*CP 50 CONTINUE C 3703 CONTINUE ENDIF DO 70 I=1,LENC HOL(I,IC) = TX2(I) QS1(I) = ALF(I,IC) + BET(I,IC)*POI(I,IC) QOL(I,IC) = MIN(QS1(I)*RHMAX,QOI(I,IC)) c TEM1 = TX2(I) + PRJ(I,IC1) * POI(I,IC) * CP HOL(I,IC) = TEM1 + QOL(I,IC) * ALHL HST(I,IC) = TEM1 + QS1(I) * ALHL C TX3(I ) = (PRJ(I,IC1) - PRH(I,IC)) * POI(I,IC) ETA(I,IC) = ETA(I,IC1) + CPBG * TX3(I) 70 CONTINUE C DO 130 I=1,LENC TX2(I) = HOL(I,K) - HST(I,IC) TX1(I) = ZERO 130 CONTINUE C C ENTRAINMENT PARAMETER C DO 160 L=IC,KM1 DO 160 I=1,LENC TX1(I) = TX1(I) + (HST(I,IC) - HOL(I,L)) * (ETA(I,L) - ETA(I,L+1)) 160 CONTINUE C LEN1 = 0 LEN2 = 0 ISAV = 0 DO 195 I=1,LENC IF (TX1(I) .GT. ZERO .AND. TX2(I) .GT. ZERO . .AND. rhfrac(i).ne.0.0 ) THEN LEN1 = LEN1 + 1 IA(LEN1) = I ALM(LEN1) = TX2(I) / TX1(I) ENDIF 195 CONTINUE C LEN2 = LEN1 if (IC1 .lt. K) then DO 196 I=1,LENC IF (TX2(I) .LE. 0.0 .AND. (HOL(I,K) .GT. HST(I,IC1)) . .AND. rhfrac(i).ne.0.0 ) THEN LEN2 = LEN2 + 1 IA(LEN2) = I ALM(LEN2) = 0.0 ENDIF 196 CONTINUE endif C IF (LEN2 .EQ. 0) THEN DO 5010 I=1,LENC*K HST(I,1) = 0.0 QOL(I,1) = 0.0 5010 CONTINUE DO 5020 I=1,LENC PCU(I) = 0.0 5020 CONTINUE RETURN ENDIF LEN11 = LEN1 + 1 C C NORMALIZED MASSFLUX C DO 250 I=1,LEN2 ETA(I,K) = 1.0 II = IA(I) TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1)) TX4(I) = PRS(II,K) 250 CONTINUE C DO 252 I=LEN11,LEN2 WFN(I) = 0.0 II = IA(I) IF (HST(II,IC1) .LT. HST(II,IC)) THEN TX6(I) = (HST(II,IC1)-HOL(II,K))/(HST(II,IC1)-HST(II,IC)) ELSE TX6(I) = 0.0 ENDIF TX2(I) = 0.5 * (PRS(II,IC1)+PRS(II,IC1+1)) * (1.0-TX6(I)) * + TX2(I) * TX6(I) 252 CONTINUE C CALL ACRITN(LEN2, TX2, TX4, TX3) C DO 260 L=KM1,IC,-1 DO 255 I=1,LEN2 TX1(I) = ETA(IA(I),L) 255 CONTINUE DO 260 I=1,LEN2 ETA(I,L) = 1.0 + ALM(I) * TX1(I) 260 CONTINUE C C CLOUD WORKFUNCTION C IF (LEN1 .GT. 0) THEN DO 270 I=1,LEN1 II = IA(I) WFN(I) = - GAM(II,IC) * (PRJ(II,IC1) - PRH(II,IC)) * * HST(II,IC) * ETA(I,IC1) 270 CONTINUE ENDIF C DO 290 I=1,LEN2 II = IA(I) TX1(I) = HOL(II,K) 290 CONTINUE C IF (IC1 .LE. KM1) THEN DO 380 L=KM1,IC1,-1 DO 380 I=1,LEN2 II = IA(I) TEM = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * HOL(II,L) C PCU(I) = PRJ(II,L+1) - PRH(II,L) TEM1 = ETA(I,L+1) * PCU(I) TX1(I) = TX1(I)*PCU(I) C PCU(I) = PRH(II,L) - PRJ(II,L) TEM1 = (TEM1 + ETA(I,L) * PCU(I)) * HST(II,L) TX1(I) = TX1(I) + TEM*PCU(I) C WFN(I) = WFN(I) + (TX1(I) - TEM1) * GAM(II,L) TX1(I) = TEM 380 CONTINUE ENDIF C LENA = 0 IF (LEN1 .GT. 0) THEN DO 512 I=1,LEN1 II = IA(I) WFN(I) = WFN(I) + TX1(I) * GAM(II,IC)*(PRJ(II,IC1)-PRH(II,IC)) * - TX3(I) IF (WFN(I) .GT. 0.0) THEN LENA = LENA + 1 I1(LENA) = IA(I) I2(LENA) = I TX1(LENA) = WFN(I) TX2(LENA) = QS1(IA(I)) TX6(LENA) = 1.0 ENDIF 512 CONTINUE ENDIF LENB = LENA DO 515 I=LEN11,LEN2 WFN(I) = WFN(I) - TX3(I) IF (WFN(I) .GT. 0.0 .AND. TX6(I) .GT. 0.0) THEN LENB = LENB + 1 I1(LENB) = IA(I) I2(LENB) = I TX1(LENB) = WFN(I) TX2(LENB) = QS1(IA(I)) TX4(LENB) = TX6(I) ENDIF 515 CONTINUE C IF (LENB .LE. 0) THEN DO 5030 I=1,LENC*K HST(I,1) = 0.0 QOL(I,1) = 0.0 5030 CONTINUE DO 5040 I=1,LENC PCU(I) = 0.0 5040 CONTINUE RETURN ENDIF C DO 516 I=1,LENB WFN(I) = TX1(I) QS1(I) = TX2(I) 516 CONTINUE C DO 520 L=IC,K DO 517 I=1,LENB TX1(I) = ETA(I2(I),L) 517 CONTINUE DO 520 I=1,LENB ETA(I,L) = TX1(I) 520 CONTINUE C LENA1 = LENA + 1 C DO 510 I=1,LENA II = I1(I) TX8(I) = HST(II,IC) - HOL(II,IC) 510 CONTINUE DO 530 I=LENA1,LENB II = I1(I) TX6(I) = TX4(I) TEM = TX6(I) * (HOL(II,IC)-HOL(II,IC1)) + HOL(II,IC1) TX8(I) = HOL(II,K) - TEM TEM1 = TX6(I) * (QOL(II,IC)-QOL(II,IC1)) + QOL(II,IC1) TX5(I) = TEM - TEM1 * ALHL QS1(I) = TEM1 + TX8(I)*(ONE/ALHL) TX3(I) = HOL(II,IC) 530 CONTINUE C C DO 620 I=1,LENB II = I1(I) WLQ(I) = QOL(II,K) - QS1(I) * ETA(I,IC) TX7(I) = HOL(II,K) 620 CONTINUE DO NT=1,Ntracer DO 621 I=1,LENB II = I1(I) UHT(I,NT) = UOI(II,K+nltop-1,NT)-UOI(II,IC+nltop-1,NT) * ETA(I,IC) 621 CONTINUE ENDDO C DO 635 L=KM1,IC,-1 DO 630 I=1,LENB II = I1(I) TEM = ETA(I,L) - ETA(I,L+1) WLQ(I) = WLQ(I) + TEM * QOL(II,L) 630 CONTINUE 635 CONTINUE DO NT=1,Ntracer DO L=KM1,IC,-1 DO I=1,LENB II = I1(I) TEM = ETA(I,L) - ETA(I,L+1) UHT(I,NT) = UHT(I,NT) + TEM * UOI(II,L+nltop-1,NT) ENDDO ENDDO ENDDO C C CALCULATE GS AND PART OF AKM (THAT REQUIRES ETA) C DO 690 I=1,LENB II = I1(I) c TX7(I) = HOL(II,K) TEM = (POI(II,KM1) - POI(II,K)) / (PRH(II,K) - PRH(II,KM1)) HOL(I,K) = TEM * (PRJ(II,K)-PRH(II,KM1))*PRH(II,K)*PRI(II,K) HOL(I,KM1) = TEM * (PRH(II,K)-PRJ(II,K))*PRH(II,KM1)*PRI(II,KM1) AKM(I) = ZERO TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1)) 690 CONTINUE IF (IC1 .LE. KM1) THEN DO 750 L=KM1,IC1,-1 DO 750 I=1,LENB II = I1(I) TEM = (POI(II,L-1) - POI(II,L)) * ETA(I,L) * / (PRH(II,L) - PRH(II,L-1)) C HOL(I,L) = TEM * (PRJ(II,L)-PRH(II,L-1)) * PRH(II,L) * * PRI(II,L) + HOL(I,L) HOL(I,L-1) = TEM * (PRH(II,L)-PRJ(II,L)) * PRH(II,L-1) * * PRI(II,L-1) C AKM(I) = AKM(I) - HOL(I,L) * * (ETA(I,L) * (PRH(II,L)-PRJ(II,L)) + * ETA(I,L+1) * (PRJ(II,L+1)-PRH(II,L))) / PRH(II,L) 750 CONTINUE ENDIF C C CALL RNCL(LENB, TX2, TX1, CLF) DO 770 I=1,LENB TX2(I) = (ONE - TX1(I)) * WLQ(I) WLQ(I) = TX1(I) * WLQ(I) C TX1(I) = HOL(I,IC) 770 CONTINUE DO 790 I=LENA1, LENB II = I1(I) TX1(I) = TX1(I) + (TX5(I)-TX3(I)+QOL(II,IC)*ALHL)*(PRI(II,IC)/CP) 790 CONTINUE DO 800 I=1,LENB HOL(I,IC) = TX1(I) - TX2(I) * ALBCP * PRI(I1(I),IC) 800 CONTINUE IF (LENA .GT. 0) THEN DO 810 I=1,LENA II = I1(I) AKM(I) = AKM(I) - ETA(I,IC1) * (PRJ(II,IC1) - PRH(II,IC)) * * TX1(I) / PRH(II,IC) 810 CONTINUE ENDIF c C CALCULATE GH C DO 830 I=1,LENB II = I1(I) TX3(I) = QOL(II,KM1) - QOL(II,K) GMH(I,K) = HOL(I,K) + TX3(I) * PRI(II,K) * (ALBCP) AKM(I) = AKM(I) + GAM(II,KM1)*(PRJ(II,K)-PRH(II,KM1)) * * GMH(I,K) TX3(I) = zero 830 CONTINUE C IF (IC1 .LE. KM1) THEN DO 840 L=KM1,IC1,-1 DO 840 I=1,LENB II = I1(I) TX2(I) = TX3(I) TX3(I) = (QOL(II,L-1) - QOL(II,L)) * ETA(I,L) TX2(I) = TX2(I) + TX3(I) C GMH(I,L) = HOL(I,L) + TX2(I) * PRI(II,L) * (ALBCP*HALF) 840 CONTINUE C C ENDIF DO 850 I=LENA1,LENB TX3(I) = TX3(I) + TWOBAL * * (TX7(I) - TX8(I) - TX5(I) - QOL(I1(I),IC)*ALHL) 850 CONTINUE DO 860 I=1,LENB GMH(I,IC) = TX1(I) + PRI(I1(I),IC) * ONEBCP * * (TX3(I)*(ALHL*HALF) + ETA(I,IC) * TX8(I)) 860 CONTINUE C C CALCULATE HC PART OF AKM C IF (IC1 .LE. KM1) THEN DO 870 I=1,LENB TX1(I) = GMH(I,K) 870 CONTINUE DO 3725 L=KM1,IC1,-1 DO 880 I=1,LENB II = I1(I) TX1(I) = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * GMH(I,L) TX2(I) = GAM(II,L-1) * (PRJ(II,L) - PRH(II,L-1)) 880 CONTINUE C IF (L .EQ. IC1) THEN DO 890 I=LENA1,LENB TX2(I) = ZERO 890 CONTINUE ENDIF DO 900 I=1,LENB II = I1(I) AKM(I) = AKM(I) + TX1(I) * * (TX2(I) + GAM(II,L)*(PRH(II,L)-PRJ(II,L))) 900 CONTINUE 3725 CONTINUE ENDIF C DO 920 I=LENA1,LENB II = I1(I) TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1)) * + 0.5*(PRS(II,IC+2) - PRS(II,IC)) * (ONE-TX6(I)) c TX1(I) = PRS(II,IC1) TX5(I) = 0.5 * (PRS(II,IC1) + PRS(II,IC+2)) C IF ((TX2(I) .GE. TX1(I)) .AND. (TX2(I) .LT. TX5(I))) THEN TX6(I) = ONE - (TX2(I) - TX1(I)) / (TX5(I) - TX1(I)) C TEM = PRI(II,IC1) / PRI(II,IC) HOL(I,IC1) = HOL(I,IC1) + HOL(I,IC) * TEM HOL(I,IC) = ZERO C GMH(I,IC1) = GMH(I,IC1) + GMH(I,IC) * TEM GMH(I,IC) = ZERO ELSEIF (TX2(I) .LT. TX1(I)) THEN TX6(I) = 1.0 ELSE TX6(I) = 0.0 ENDIF 920 CONTINUE C C DO I=1,LENC PCU(I) = 0.0 ENDDO DO 970 I=1,LENB II = I1(I) IF (AKM(I) .LT. ZERO .AND. WLQ(I) .GE. 0.0) THEN WFN(I) = - TX6(I) * WFN(I) * RASALF / AKM(I) ELSE WFN(I) = ZERO ENDIF TEM = (PRS(II,K+1)-PRS(II,K))*(CMB2PA*FRAC) WFN(I) = MIN(WFN(I), TEM) C C compute cloud amount C CC TX1(I) = CLN(II) CC IF (WFN(I) .GT. CRTMSF) TX1(I) = TX1(I) + CLF(I) CC IF (TX1(I) .GT. ONE) TX1(I) = ONE C C PRECIPITATION C PCU(II) = WLQ(I) * WFN(I) * ONEBG C C CUMULUS FRICTION AT THE BOTTOM LAYER C TX4(I) = WFN(I) * (1.0/ALHL) TX5(I) = WFN(I) * ONEBCP 970 CONTINUE C C compute cloud mass flux for diagnostic output C DO L = IC,K DO I=1,LENB II = I1(I) if(L.lt.K)then CMASS(II,L) = ETA(I,L+1) * WFN(I) * ONEBG else CMASS(II,L) = WFN(I) * ONEBG endif ENDDO ENDDO C CC DO 975 I=1,LENB CC II = I1(I) CC CLN(II) = TX1(I) CC975 CONTINUE C C THETA AND Q CHANGE DUE TO CLOUD TYPE IC C c TEMA = 0.0 c TEMB = 0.0 DO 990 L=IC,K DO 980 I=1,LENB II = I1(I) TEM = (GMH(I,L) - HOL(I,L)) * TX4(I) TEM1 = HOL(I,L) * TX5(I) C TCU(II,L) = TEM1 / PRH(II,L) QCU(II,L) = TEM 980 CONTINUE c I = I1(IP1) c c TEM = (PRS(I,L+1)-PRS(I,L)) * (ONEBG*100.0) c TEMA = TEMA + TCU(I,L) * PRH(I,L) * TEM * (CP/ALHL) c TEMB = TEMB + QCU(I,L) * TEM C 990 CONTINUE C c Compute Tracer Tendencies c ------------------------- do nt = 1,ntracer c Tracer Tendency at the Bottom Layer c ----------------------------------- DO 995 I=1,LENB II = I1(I) TEM = half*TX5(I) * PRI(II,K) TX1(I) = (UOI(II,KM1+nltop-1,nt) - UOI(II,K+nltop-1,nt)) ucu(II,K,nt) = TEM * TX1(I) 995 CONTINUE c Tracer Tendency at all other Levels c ----------------------------------- DO 1020 L=KM1,IC1,-1 DO 1010 I=1,LENB II = I1(I) TEM = half*TX5(I) * PRI(II,L) TEM1 = TX1(I) TX1(I) = (UOI(II,L-1+nltop-1,nt)-UOI(II,L+nltop-1,nt)) * ETA(I,L) TX3(I) = (TX1(I) + TEM1) * TEM 1010 CONTINUE DO 1020 I=1,LENB II = I1(I) ucu(II,L,nt) = TX3(I) 1020 CONTINUE DO 1030 I=1,LENB II = I1(I) IF (TX6(I) .GE. 1.0) THEN TEM = half*TX5(I) * PRI(II,IC) ELSE TEM = 0.0 ENDIF TX1(I) = (TX1(I) + UHT(I,nt) + UHT(I,nt)) * TEM 1030 CONTINUE DO 1040 I=1,LENB II = I1(I) ucu(II,IC,nt) = TX1(I) 1040 CONTINUE enddo C C PENETRATIVE CONVECTION CALCULATION OVER C RETURN END SUBROUTINE RNCL(LEN, PL, RNO, CLF) C C C********************************************************************* C********************** Relaxed Arakawa-Schubert ********************* C************************ SUBROUTINE RNCL ************************ C**************************** 23 July 1992 *************************** C********************************************************************* PARAMETER (P5=500.0, P8=800.0, PT8=0.8, PT2=0.2) PARAMETER (PFAC=PT2/(P8-P5)) C PARAMETER (P4=400.0, P6=401.0) PARAMETER (P7=700.0, P9=900.0) PARAMETER (CUCLD=0.5,CFAC=CUCLD/(P6-P4)) C DIMENSION PL(LEN), RNO(LEN), CLF(LEN) DO 10 I=1,LEN rno(i) = 1.0 ccc if( pl(i).le.400.0 ) rno(i) = max( 0.75, 1.0-0.0025*(400.0-pl(i)) ) ccc IF ( PL(I).GE.P7 .AND. PL(I).LE.P9 ) THEN ccc RNO(I) = ((P9-PL(I))/(P9-P7)) **2 ccc ELSE IF (PL(I).GT.P9) THEN ccc RNO(I) = 0. ccc ENDIF CLF(I) = CUCLD C CARIESIF (PL(I) .GE. P5 .AND. PL(I) .LE. P8) THEN CARIES RNO(I) = (P8-PL(I))*PFAC + PT8 CARIESELSEIF (PL(I) .GT. P8 ) THEN CARIES RNO(I) = PT8 CARIESENDIF CARIES IF (PL(I) .GE. P4 .AND. PL(I) .LE. P6) THEN CLF(I) = (P6-PL(I))*CFAC ELSEIF (PL(I) .GT. P6 ) THEN CLF(I) = 0.0 ENDIF 10 CONTINUE C RETURN END SUBROUTINE ACRITN ( LEN,PL,PLB,ACR ) C********************************************************************* C********************** Relaxed Arakawa-Schubert ********************* C************************** SUBROUTINE ACRIT ********************* C****************** modified August 28, 1996 L.Takacs ************ C**** ***** C**** Note: Data obtained from January Mean After-Analysis ***** C**** from 4x5 46-layer GEOS Assimilation ***** C**** ***** C********************************************************************* real PL(LEN), PLB(LEN), ACR(LEN) parameter (lma=18) real p(lma) real a(lma) data p / 93.81, 111.65, 133.46, 157.80, 186.51, . 219.88, 257.40, 301.21, 352.49, 409.76, . 471.59, 535.04, 603.33, 672.79, 741.12, . 812.52, 875.31, 930.20/ data a / 3.35848, 3.13645, 2.48072, 2.08277, 1.53364, . 1.01971, .65846, .45867, .38687, .31002, . .25574, .20347, .17254, .15260, .16756, . .09916, .10360, .05880/ do L=1,lma-1 do i=1,len if( pl(i).ge.p(L) .and. . pl(i).le.p(L+1)) then temp = ( pl(i)-p(L) )/( p(L+1)-p(L) ) acr(i) = a(L+1)*temp + a(L)*(1-temp) endif enddo enddo do i=1,len if( pl(i).lt.p(1) ) acr(i) = a(1) if( pl(i).gt.p(lma) ) acr(i) = a(lma) enddo do i=1,len acr(i) = acr(i) * (plb(i)-pl(i)) enddo RETURN END SUBROUTINE RNEVP(NN,IRUN,NLAY,TL,QL,RAIN,PL,CLFRAC,SP,DSIG,PLKE, 1 PLK,TH,TEMP1,TEMP2,TEMP3,ITMP1,ITMP2,RCON,RLAR,CLSBTH,tmscl, 2 tmfrc,cp,gravity,alhl,gamfac,cldlz,RHCRIT,offset,alpha) PARAMETER (ZM1P04 = -1.04E-4 ) PARAMETER (ZERO = 0.) PARAMETER (TWO89= 2.89E-5) PARAMETER ( ZP44= 0.44) PARAMETER ( ZP01= 0.01) PARAMETER ( ZP1 = 0.1 ) PARAMETER ( ZP001= 0.001) PARAMETER ( HALF= 0.5) PARAMETER ( ZP578 = 0.578 ) PARAMETER ( ONE = 1.) PARAMETER ( THOUSAND = 1000.) PARAMETER ( Z3600 = 3600.) C DIMENSION TL(IRUN,NLAY),QL(IRUN,NLAY),RAIN(IRUN,NLAY), $ PL(IRUN,NLAY),CLFRAC(IRUN,NLAY),SP(IRUN),TEMP1(IRUN,NLAY), $ TEMP2(IRUN,NLAY),PLKE(IRUN,NLAY), $ RCON(IRUN),RLAR(IRUN),DSIG(NLAY),PLK(IRUN,NLAY),TH(IRUN,NLAY), $ TEMP3(IRUN,NLAY),ITMP1(IRUN,NLAY), $ ITMP2(IRUN,NLAY),CLSBTH(IRUN,NLAY) C DIMENSION EVP9(IRUN,NLAY) real water(irun),crystal(irun) real watevap(irun),iceevap(irun) real fracwat,fracice, tice,rh,fact,dum real cldlz(irun,nlay) real rhcrit(irun,nlay), rainmax(irun) real offset, alpha c Explicit Inline Directives c -------------------------- #if CRAY #if f77 cfpp$ expand (qsat) #endif #endif tice = getcon('FREEZING-POINT') fracwat = 0.70 fracice = 0.01 NLAYM1 = NLAY - 1 IRNLAY = IRUN*NLAY IRNLM1 = IRUN*(NLAY-1) RPHF = Z3600/tmscl ELOCP = alhl/cp CPOG = cp/gravity DO I = 1,IRUN RLAR(I) = 0. water(i) = 0. crystal(i) = 0. ENDDO do L = 1,nlay do i = 1,irun EVP9(i,L) = 0. TEMP1(i,L) = 0. TEMP2(i,L) = 0. TEMP3(i,L) = 0. CLSBTH(i,L) = 0. cldlz(i,L) = 0. enddo enddo C RHO(ZERO) / RHO FOR TERMINAL VELOCITY APPROX. c --------------------------------------------- DO L = 1,NLAY DO I = 1,IRUN TEMP2(I,L) = PL(I,L)*ZP001 TEMP2(I,L) = SQRT(TEMP2(I,L)) ENDDO ENDDO C INVERSE OF MASS IN EACH LAYER c ----------------------------- DO L = 1,NLAY DO I = 1,IRUN TEMP3(I,L) = SP(I) * DSIG(L) TEMP3(I,L) = GRAVITY*ZP01 / TEMP3(I,L) ENDDO ENDDO C DO LOOP FOR MOISTURE EVAPORATION ABILITY AND CONVEC EVAPORATION. c ---------------------------------------------------------------- DO 100 L=1,NLAY DO I = 1,IRUN TEMP1(I,3) = TL(I,L) TEMP1(I,4) = QL(I,L) ENDDO DO 50 N=1,2 IF(N.EQ.1)RELAX=HALF IF(N.GT.1)RELAX=ONE DO I = 1,IRUN call qsat ( temp1(i,3),pl(i,L),temp1(i,2),temp1(i,6),.true. ) TEMP1(I,5)=TEMP1(I,2)-TEMP1(I,4) TEMP1(I,6)=TEMP1(I,6)*ELOCP TEMP1(I,5)=TEMP1(I,5)/(ONE+TEMP1(I,6)) TEMP1(I,4)=TEMP1(I,4)+TEMP1(I,5)*RELAX TEMP1(I,3)=TEMP1(I,3)-TEMP1(I,5)*ELOCP*RELAX ENDDO 50 CONTINUE DO I = 1,IRUN EVP9(I,L) = (TEMP1(I,4) - QL(I,L))/TEMP3(I,L) C convective detrained water cldlz(i,L) = rain(i,L)*temp3(i,L) if( tl(i,L).gt.tice-20.) then water(i) = water(i) + rain(i,L) else crystal(i) = crystal(i) + rain(i,L) endif ENDDO C********************************************************************** C FOR CONVECTIVE PRECIP, FIND THE "EVAPORATION EFFICIENCY" USING * C KESSLERS PARAMETERIZATION * C********************************************************************** DO 20 I=1,IRUN iceevap(i) = 0. watevap(i) = 0. if( (evp9(i,L).gt.0.) .and. (crystal(i).gt.0.) ) then iceevap(I) = EVP9(I,L)*fracice IF(iceevap(i).GE.crystal(i)) iceevap(i) = crystal(i) EVP9(I,L)=EVP9(I,L)-iceevap(I) crystal(i) = crystal(i) - iceevap(i) endif C and now warm precipitate if( (evp9(i,L).gt.0.) .and. (water(i).gt.0.) ) then exparg = ZM1P04*tmscl*((water(i)*RPHF*TEMP2(I,L))**ZP578) AREARAT = ONE-(EXP(EXPARG)) watevap(I) = EVP9(I,L)*AREARAT*fracwat IF(watevap(I).GE.water(i)) watevap(I) = water(i) EVP9(I,L)=EVP9(I,L)-watevap(I) water(i) = water(i) - watevap(i) endif QL(I,L) = QL(I,L)+(iceevap(i)+watevap(i))*TEMP3(I,L) TL(I,L) = TL(I,L)-(iceevap(i)+watevap(i))*TEMP3(I,L)*ELOCP 20 CONTINUE 100 CONTINUE do i = 1,irun rcon(i) = water(i) + crystal(i) enddo C********************************************************************** C Large Scale Precip C********************************************************************** DO 200 L=1,NLAY DO I = 1,IRUN rainmax(i) = rhcrit(i,L)*evp9(i,L) + . ql(i,L)*(rhcrit(i,L)-1.)/temp3(i,L) if (rainmax(i).LE.0.0) then call qsat( tl(i,L),pl(i,L),rh,dum,.false.) rh = ql(i,L)/rh if( rhcrit(i,L).eq.1.0 ) then fact = 1.0 else fact = min( 1.0, alpha + (1.0-alpha)*( rh-rhcrit(i,L)) / 1 (1.0-rhcrit(i,L)) ) endif C Do not allow clouds above 10 mb if( pl(i,L).ge.10.0 ) CLSBTH(I,L) = fact RLAR(I) = RLAR(I)-rainmax(I) QL(I,L) = QL(I,L)+rainmax(I)*TEMP3(I,L) TL(I,L) = TL(I,L)-rainmax(I)*TEMP3(I,L)*ELOCP C Large-scale water cldlz(i,L) = cldlz(i,L) - rainmax(i)*temp3(i,L) ENDIF ENDDO DO I=1,IRUN IF((RLAR(I).GT.0.0).AND.(rainmax(I).GT.0.0))THEN RPOW=(RLAR(I)*RPHF*TEMP2(I,L))**ZP578 EXPARG = ZM1P04*tmscl*RPOW AREARAT = ONE-(EXP(EXPARG)) TEMP1(I,7) = rainmax(I)*AREARAT IF(TEMP1(I,7).GE.RLAR(I)) TEMP1(I,7) = RLAR(I) RLAR(I) = RLAR(I)-TEMP1(I,7) QL(I,L) = QL(I,L)+TEMP1(I,7)*TEMP3(I,L) TL(I,L) = TL(I,L)-TEMP1(I,7)*TEMP3(I,L)*ELOCP ENDIF ENDDO 200 CONTINUE RETURN END subroutine srclouds (th,q,plk,pl,plke,cloud,cldwat,irun,irise, 1 rhc,offset,alpha) C*********************************************************************** C C PURPOSE: C ======== C Compute non-precipitating cloud fractions C based on Slingo and Ritter (1985). C Remove cloudiness where conditionally unstable. C C INPUT: C ====== C th ......... Potential Temperature (irun,irise) C q .......... Specific Humidity (irun,irise) C plk ........ P**Kappa at mid-layer (irun,irise) C pl ......... Pressure at mid-layer (irun,irise) C plke ....... P**Kappa at edge (irun,irise+1) C irun ....... Horizontal dimension C irise ...... Vertical dimension C C OUTPUT: C ======= C cloud ...... Cloud Fraction (irun,irise) C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** implicit none integer irun,irise real th(irun,irise) real q(irun,irise) real plk(irun,irise) real pl(irun,irise) real plke(irun,irise+1) real tempth(irun) real tempqs(irun) real dhstar(irun) real cloud(irun,irise) real cldwat(irun,irise) real qs(irun,irise) real cp, alhl, getcon, akap, pcheck real ratio, temp, pke, elocp real rhcrit,rh,dum,pbar,tbar integer i,L,ntradesu,ntradesl real factor real rhc(irun,irise) real offset,alpha c Explicit Inline Directives c -------------------------- #if CRAY #if f77 cfpp$ expand (qsat) #endif #endif cp = getcon('CP') alhl = getcon('LATENT HEAT COND') elocp = alhl/cp akap = getcon('KAPPA') do L = 1,irise do i = 1,irun temp = th(i,L)*plk(i,L) call qsat ( temp,pl(i,L),qs(i,L),dum,.false. ) enddo enddo do L = 2,irise do i = 1,irun rh = q(i,L)/qs(i,L) rhcrit = rhc(i,L) - offset ratio = alpha*(rh-rhcrit)/offset if(cloud(i,L).eq. 0.0 .and. ratio.gt.0.0 ) then cloud(i,L) = min( ratio,1.0 ) endif enddo enddo c Reduce clouds from conditionally unstable layer c ----------------------------------------------- call ctei ( th,q,cloud,cldwat,pl,plk,plke,irun,irise ) return end subroutine ctei ( th,q,cldfrc,cldwat,pl,plk,plke,im,lm ) implicit none integer im,lm real th(im,lm),q(im,lm),plke(im,lm+1),cldwat(im,lm) real plk(im,lm),pl(im,lm),cldfrc(im,lm) integer i,L real getcon,cp,alhl,elocp,cpoel,t,p,s,qs,dqsdt,dq real k,krd,kmm,f cp = getcon('CP') alhl = getcon('LATENT HEAT COND') cpoel = cp/alhl elocp = alhl/cp do L=lm,2,-1 do i=1,im dq = q(i,L)+cldwat(i,L)-q(i,L-1)-cldwat(i,L-1) if( dq.eq.0.0 ) dq = 1.0e-20 k = 1.0 + cpoel*plke(i,L)*( th(i,L)-th(i,L-1) ) / dq t = th(i,L)*plk(i,L) p = pl(i,L) call qsat ( t,p,qs,dqsdt,.true. ) krd = ( cpoel*t*(1+elocp*dqsdt) )/( 1 + 1.608*dqsdt*t ) kmm = ( 1+elocp*dqsdt )*( 1 + 0.392*cpoel*t ) . / ( 2+(1+1.608*cpoel*t)*elocp*dqsdt ) s = ( (k-krd)/(kmm-krd) ) f = 1.0 - min( 1.0, max(0.0,1.0-exp(-s)) ) cldfrc(i,L) = cldfrc(i,L)*f cldwat(i,L) = cldwat(i,L)*f enddo enddo return end subroutine back2grd(gathered,indeces,scattered,irun) implicit none integer i,irun,indeces(irun) real gathered(irun),scattered(irun) real temp(irun) do i = 1,irun temp(indeces(i)) = gathered(i) enddo do i = 1,irun scattered(i) = temp(i) enddo return end