/[MITgcm]/MITgcm/pkg/fizhi/fizhi_turb.F
ViewVC logotype

Diff of /MITgcm/pkg/fizhi/fizhi_turb.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by molod, Tue Jun 15 16:06:03 2004 UTC revision 1.23 by molod, Tue Aug 10 15:13:47 2004 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "FIZHI_OPTIONS.h"
5        subroutine turbio (im,jm,nlay,maxtypes,istrip,nymd,nhms        subroutine turbio (im,jm,nlay,istrip,nymd,nhms,bi,bj
6       1 ,ndturb,ptop, pz, uz, vz, tz, qz, ntracers,ptracers,       1 ,ndturb,ptop, pz, uz, vz, tz, qz, ntracers,ptracers
7       2 ,plz,plze,dpres,pkht,pkz,ctmt,xxmt,yymt,zetamt,xlmt,khmt,tke       2 ,plz,plze,dpres,pkht,pkz,ctmt,xxmt,yymt,zetamt,xlmt,khmt,tke
8       3 ,tgz,gwet,snow,fracland,landtype       3 ,tgz,fracland,landtype
9       4 ,tcanopy,ecanopy,tdeep,swetshal,swetroot,swetdeep,snodep,capac       4 ,tcanopy,ecanopy,tdeep,swetshal,swetroot,swetdeep,snodep,capac
10       5 ,nchp,nchplnd,chfr,chlt,chlon,igrd,ityp,alai,agrn       5 ,nchp,nchptot,nchplnd,chfr,chlt,chlon,igrd,ityp,alai,agrn,thkz
11       6 ,surftype,tilefrac,thkz,tprof       6 ,tprof
12       8 ,duturb, dvturb, dtturb,dqturb,radlwg,st4,dst4,radswg,radswt       8 ,duturb, dvturb, dtturb,dqturb,radlwg,st4,dst4,radswg,radswt
13       9 ,fdifpar,fdirpar,rainlsp,rainconv,snowfall,tempref       9 ,fdifpar,fdirpar,rainlsp,rainconv,snowfall,tempref
14       1 ,imstturblw,imstturbsw,qliqavelw,qliqavesw,fccavelw,fccavesw       1 ,imstturblw,imstturbsw,qliqavelw,qliqavesw,fccavelw,fccavesw
15       2 ,qqgrid,u2m,v2m,t2m,q2m,u10m,v10m,t10m,q10m,myid,comm)       2 ,qqgrid,myid)
16  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
17  c subroutine turbio - model interface routine to trbflx, the turbulence  c subroutine turbio - model interface routine to trbflx, the turbulence
18  c        parameterization, and tile, the land surface parameterization  c        parameterization, and tile, the land surface parameterization
# Line 21  c  input: Line 21  c  input:
21  c      im      - number of points in the longitude direction  c      im      - number of points in the longitude direction
22  c      jm      - number of points in the latitude direction  c      jm      - number of points in the latitude direction
23  c      nlay    - number of vertical levels  c      nlay    - number of vertical levels
 c      maxtypes- maximum number of different surface types in any one grid cell  
24  c      istrip  - number of horizontal points to be handled at a time on  c      istrip  - number of horizontal points to be handled at a time on
25  c      nymd    - year and date integer in YYMMDD format (ie, 790212)  c      nymd    - year and date integer in YYMMDD format (ie, 790212)
26  c      nhms    - date and time integer in HHMMSS format (ie, 123000)  c      nhms    - date and time integer in HHMMSS format (ie, 123000)
# Line 46  c      igrd    - tile space grid number Line 45  c      igrd    - tile space grid number
45  c      ityp    - tile space vegetation type - integer[nchp]  c      ityp    - tile space vegetation type - integer[nchp]
46  c      alai    - leaf area index - real[nchp]  c      alai    - leaf area index - real[nchp]
47  c      agrn    - greenness fraction - real[nchp]  c      agrn    - greenness fraction - real[nchp]
 c      surftype- surface type designation - integer[lon,lat,maxtype]  
 c                where maxtype is the maximum number of tiles in a grid cell  
 c      tilefrac- fraction of grid box covered by the corresponding surface type  
 c                designation - real[lon,lat,maxtype]  
48  c      thkz    - sea ice thickness in m (0. for no ice) - real[lon,lat]  c      thkz    - sea ice thickness in m (0. for no ice) - real[lon,lat]
49  c      tprof   - logical flag for point by point diagnostic output  c      tprof   - logical flag for point by point diagnostic output
50  c      ndiagsiz- number of diagnostic 2-D arrays allocated  c      ndiagsiz- number of diagnostic 2-D arrays allocated
# Line 70  c      snowfall- total snowfall rate in Line 65  c      snowfall- total snowfall rate in
65  c updated:  c updated:
66  c      tke     - turbulent k.e. in m**2/s**2 - real[tiles,levels]  c      tke     - turbulent k.e. in m**2/s**2 - real[tiles,levels]
67  c      tgz     - surface skin temperature in deg K - real[lon,lat]  c      tgz     - surface skin temperature in deg K - real[lon,lat]
 c      gwet    - surface ground wetness fraction - real[lon,lat]  
 c      snow    - depth of snow pack in cm liquid water equiv real[lon,lat]  
68  c      tcanopy - canopy temperature in deg K real[tiles]  c      tcanopy - canopy temperature in deg K real[tiles]
69  c                 (sea surface temp over the ocean tiles)  c                 (sea surface temp over the ocean tiles)
70  c      ecanopy - canopy vapor pressure in mb real[tiles]  c      ecanopy - canopy vapor pressure in mb real[tiles]
# Line 96  c      qliqavesw - Moist   Turbulence Li Line 89  c      qliqavesw - Moist   Turbulence Li
89  c      fccavelw  - Moist   Turbulence Cloud Fraction for Longwave  - real[lon,lat,levels]  c      fccavelw  - Moist   Turbulence Cloud Fraction for Longwave  - real[lon,lat,levels]
90  c      fccavesw  - Moist   Turbulence Cloud Fraction for Shortwave - real[lon,lat,levels]  c      fccavesw  - Moist   Turbulence Cloud Fraction for Shortwave - real[lon,lat,levels]
91  c      qqgrid    - Gridded Turbulent Kinetic Energy                - real[lon,lat,levels]  c      qqgrid    - Gridded Turbulent Kinetic Energy                - real[lon,lat,levels]
 c      u2m       - Gridded U-Wind      real[lon,lat] at  2 meters  
 c      v2m       - Gridded V-Wind      real[lon,lat] at  2 meters  
 c      t2m       - Gridded Temperature real[lon,lat] at  2 meters  
 c      q2m       - Gridded Spec.Humd.  real[lon,lat] at  2 meters  
 c      u10m      - Gridded U-Wind      real[lon,lat] at 10 meters  
 c      v10m      - Gridded V-Wind      real[lon,lat] at 10 meters  
 c      t10m      - Gridded Temperature real[lon,lat] at 10 meters  
 c      q10m      - Gridded Spec.Humd.  real[lon,lat] at 10 meters  
 c  
92  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
93        implicit none        implicit none
94    
95    #ifdef ALLOW_USE_MPI
96  #include "mpif.h"  #include "mpif.h"
97    #endif
98    
99  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
100    #include "SIZE.h"
101    #include "diagnostics_SIZE.h"
102  #include "diagnostics.h"  #include "diagnostics.h"
103  #endif  #endif
104    
105        integer n,im,jm,nlay,maxtypes,istrip,nymd,nhms,ndturb        integer im,jm,nlay,istrip,nymd,nhms,bi,bj,ndturb
       real   ptop  
       integer nchp,nchplnd,landtype(im,jm),igrd(nchp)  
       integer ityp(nchp)  
       real fracland(im,jm),chfr(nchp),chlt(nchp)  
       real chlon(nchp)  
       real alai(nchp),agrn(nchp)  
         
       real pz(im,jm)  
       real uz(im,jm,nlay)  
       real vz(im,jm,nlay)  
       real tz(im,jm,nlay)  
106        integer ntracers, ptracers        integer ntracers, ptracers
107        real qz(im,jm,nlay,ntracers)        integer nchp,nchptot,nchplnd
108        real pkht(im,jm,nlay)        _RL ptop
109        real tke(nchp,nlay)        _RL pz(im,jm),uz(im,jm,nlay),vz(im,jm,nlay),tz(im,jm,nlay)
110        integer imstturblw, imstturbsw        _RL qz(im,jm,nlay,ntracers)
111        real ctmt(nchp),xxmt(nchp),yymt(nchp),zetamt(nchp)        _RL plz(im,jm,nlay),plze(im,jm,nlay+1),dpres(im,jm,nlay)
112        real xlmt(nchp,nlay),khmt(nchp,nlay)        _RL pkht(im,jm,nlay+1),pkz(im,jm,nlay)
113          _RL ctmt(nchp),xxmt(nchp),yymt(nchp),zetamt(nchp)
114        real tgz(im, jm),gwet(im,jm),snow(im,jm)        _RL xlmt(nchp,nlay),khmt(nchp,nlay),tke(nchp,nlay)
115        real  tcanopy(nchp)        _RL tgz(im, jm),fracland(im,jm)
116        real    tdeep(nchp)        integer landtype(im,jm)
117        real  ecanopy(nchp)        _RL tcanopy(nchp),tdeep(nchp),ecanopy(nchp),swetshal(nchp)
118        real swetshal(nchp)        _RL swetroot(nchp),swetdeep(nchp),snodep(nchp),capac(nchp)
119        real swetroot(nchp)        _RL chfr(nchp),chlt(nchp),chlon(nchp)
120        real swetdeep(nchp)        integer igrd(nchp),ityp(nchp)
121        real   snodep(nchp)        _RL alai(nchp),agrn(nchp),thkz(im,jm)
       real    capac(nchp)  
   
       integer surftype(im,jm,maxtypes)  
       real    tilefrac(im,jm,maxtypes)  
       real thkz(im,jm)  
   
122        logical tprof        logical tprof
123          _RL duturb(im,jm,nlay),dvturb(im,jm,nlay)
124          _RL dtturb(im,jm,nlay),dqturb(im,jm,nlay,ntracers)
125          _RL st4(im,jm),dst4(im,jm)
126          _RL radswg(im,jm),radswt(im,jm),radlwg(im,jm)
127          _RL fdifpar(im,jm),fdirpar(im,jm)
128          _RL rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)
129          _RL tempref (im,jm)
130          integer imstturblw, imstturbsw
131          _RL qliqavesw(im,jm,nlay),qliqavelw(im,jm,nlay)
132          _RL fccavelw (im,jm,nlay),fccavesw (im,jm,nlay)
133          _RL qqgrid   (im,jm,nlay)
134          integer myid
135    
136        real duturb(im,jm,nlay),dvturb(im,jm,nlay)  C Local Variables
       real dtturb(im,jm,nlay),dqturb(im,jm,nlay,ntracers)  
   
       real st4(im,jm),dst4(im,jm)  
       real radswg(im,jm),radswt(im,jm),radlwg(im,jm)  
       real fdifpar(im,jm),fdirpar(im,jm)  
       real rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)  
       real tempref (im,jm)  
   
       real qliqavesw(im,jm,nlay),qliqavelw(im,jm,nlay)  
       real fccavelw (im,jm,nlay),fccavesw (im,jm,nlay)  
       real qqgrid   (im,jm,nlay)  
   
 C Couplings for Output Purposes, Not Reproduceable  
       real u2m(im,jm), u10m(im,jm)    
       real v2m(im,jm), v10m(im,jm)    
       real t2m(im,jm), t10m(im,jm)    
       real q2m(im,jm), q10m(im,jm)    
137    
138        integer numstrips        integer numstrips
139        integer ijall        integer ijall
140        real fmu,hice,tref,pref,cti,ed        _RL fmu,hice,tref,pref,cti,ed
141  C Set fmu and ed to zero for no background diffusion  C Set fmu and ed to zero for no background diffusion
142        parameter ( fmu    = 0.00000    )        parameter ( fmu    = 0.00000    )
143        parameter ( hice   =   300.     )        parameter ( hice   =   300.     )
# Line 179  C Set fmu and ed to zero for no backgrou Line 146  C Set fmu and ed to zero for no backgrou
146        parameter ( cti    =   0.0052   )        parameter ( cti    =   0.0052   )
147        parameter ( ed     =   0.0      )        parameter ( ed     =   0.0      )
148    
149        real     pkz(im,jm,nlay)        _RL qliqtmp(im,jm,nlay)
150        real     plz(im,jm,nlay)        _RL  fcctmp(im,jm,nlay)
151        real qliqtmp(im,jm,nlay)        _RL tmpdiag(im,jm)
152        real  fcctmp(im,jm,nlay)        _RL   thtgz(im*jm)
153        real tmpdiag(im,jm)        _RL tempor1(im,jm,nlay),tempor2(im,jm,nlay),tempor3(im,jm,nlay)
154        real   thtgz(im*jm)        _RL tempor4(im,jm,nlay)
155          _RL tempsfc1(im,jm),tempsfc2(im,jm),tempsfc3(im,jm)
156    
157        integer nland        integer nland
158        real alwcoeff(nchp),blwcoeff(nchp)        _RL alwcoeff(nchp),blwcoeff(nchp)
159        real netsw(nchp)        _RL netsw(nchp)
160        real cnvprec(nchp),lsprec(nchp)        _RL cnvprec(nchp),lsprec(nchp)
161        real snowprec(nchp)        _RL snowprec(nchp)
162        real pardiff(nchp),pardirct(nchp)        _RL pardiff(nchp),pardirct(nchp)
163        real pmsc(nchp),radmsc(nchp)        _RL pmsc(nchp)
164        real netlw(nchp)        _RL netlw(nchp)
165        real sqscat(nchp), rsoil1(nchp)        _RL sqscat(nchp), rsoil1(nchp)
166        real rsoil2(nchp)        _RL rsoil2(nchp)
167        real rdc(nchp),u2fac(nchp)        _RL rdc(nchp),u2fac(nchp)
168        real z2ch(nchp)        _RL z2ch(nchp)
169        real zoch(nchp),cdrc(nchp)        _RL zoch(nchp),cdrc(nchp)
170        real cdsc(nchp)        _RL cdsc(nchp)
171        real dqsdt(nchp)        _RL dqsdt(nchp)
172        real tground(nchp),qground(nchp)        _RL tground(nchp),qground(nchp)
173        real utility(nchp)        _RL utility(nchp)
174        real    qice(nchp)        _RL    qice(nchp)
175        real   dqice(nchp)        _RL   dqice(nchp)
176    
177        real dumsc(nchp,nlay),dvmsc(nchp,nlay)        _RL dumsc(nchp,nlay),dvmsc(nchp,nlay)
178        real dtmsc(nchp,nlay),dqmsc(nchp,nlay,ntracers)        _RL dtmsc(nchp,nlay),dqmsc(nchp,nlay,ntracers)
179    
180        real shg(nchp),z0(nchp),icethk(nchp)        _RL shg(nchp),z0(nchp),icethk(nchp)
181        integer water(nchp)        integer water(nchp)
182    
183        real lats(istrip),lons(istrip),cosz(istrip),icest(istrip)        _RL lats(istrip),lons(istrip),cosz(istrip),icest(istrip)
184        real rainls(istrip),raincon(istrip),newsnow(istrip)        _RL rainls(istrip),raincon(istrip),newsnow(istrip)
185        real lwnet(istrip),pardf(istrip),pardr(istrip),swnet(istrip)        _RL pardf(istrip),pardr(istrip),swnet(istrip)
186        real hlwdwn(istrip),alwrad(istrip),blwrad(istrip),tmpnlay(istrip)        _RL hlwdwn(istrip),alwrad(istrip),blwrad(istrip)
187        real laistrip(istrip),grnstrip(istrip),z2str(istrip),cd(istrip)        _RL tmpnlay(istrip)
188        real scatstr(istrip), rs1str(istrip), rs2str(istrip)        _RL laistrip(istrip),grnstrip(istrip),z2str(istrip),cd(istrip)
189        real rdcstr(istrip),u2fstr(istrip),dqsdtstr(istrip)        _RL scatstr(istrip), rs1str(istrip), rs2str(istrip)
190        real eturb(istrip),dedqa(istrip),dedtc(istrip)        _RL rdcstr(istrip),u2fstr(istrip),dqsdtstr(istrip)
191        real hsturb(istrip),dhsdqa(istrip),dhsdtc(istrip)        _RL eturb(istrip),dedqa(istrip),dedtc(istrip)
192        real savetc(istrip),saveqa(istrip),lwstrip(istrip)        _RL hsturb(istrip),dhsdqa(istrip),dhsdtc(istrip)
193        real chfrstr(istrip),psurf(istrip),shgstr(istrip)        _RL savetc(istrip),saveqa(istrip),lwstrip(istrip)
194          _RL chfrstr(istrip),psurf(istrip),shgstr(istrip)
195        integer types(istrip),igrdstr(istrip)        integer types(istrip),igrdstr(istrip)
196        real evap(istrip),shflux(istrip),runoff(istrip),bomb(istrip)        _RL evap(istrip),shflux(istrip),runoff(istrip),bomb(istrip)
197        real eint(istrip),esoi(istrip),eveg(istrip),esno(istrip)        _RL eint(istrip),esoi(istrip),eveg(istrip),esno(istrip)
198        real smelt(istrip),hlatn(istrip),hlwup(istrip),gdrain(istrip)        _RL smelt(istrip),hlatn(istrip),hlwup(istrip),gdrain(istrip)
199        real runsrf(istrip),fwsoil(istrip),evpot(istrip)        _RL runsrf(istrip),fwsoil(istrip),evpot(istrip)
200        real strdg1(istrip),strdg2(istrip),strdg3(istrip),strdg4(istrip)        _RL strdg1(istrip),strdg2(istrip),strdg3(istrip),strdg4(istrip)
201        real strdg5(istrip),strdg6(istrip),strdg7(istrip),strdg8(istrip)        _RL strdg5(istrip),strdg6(istrip),strdg7(istrip),strdg8(istrip)
202        real strdg9(istrip),tmpstrip(istrip),qicestr(istrip)        _RL strdg9(istrip),tmpstrip(istrip),qicestr(istrip)
203        real dqicestr(istrip)        _RL dqicestr(istrip)
204    
205        real u(istrip,nlay+1), v(istrip,nlay+1), th(istrip,nlay+1)        _RL u(istrip,nlay+1), v(istrip,nlay+1), th(istrip,nlay+1)
206        real sh(istrip,nlay+1), thv(istrip,nlay+1), pe(istrip,nlay+1)        _RL sh(istrip,nlay+1), thv(istrip,nlay+1), pe(istrip,nlay+1)
207        real tracers(istrip,nlay+1,ntracers)        _RL tracers(istrip,nlay+1,ntracers)
208        real pke(istrip,nlay+1)        _RL dpstr(istrip,nlay),pke(istrip,nlay+1)
209        real pk(istrip,nlay), qq(istrip,nlay),   p(istrip,nlay)        _RL pk(istrip,nlay), qq(istrip,nlay),   p(istrip,nlay)
210        real sri(istrip,nlay), skh(istrip,nlay), skm(istrip,nlay)        _RL sri(istrip,nlay), skh(istrip,nlay), skm(istrip,nlay)
211        real stuflux(istrip,nlay), stvflux(istrip,nlay)        _RL stuflux(istrip,nlay), stvflux(istrip,nlay)
212        real sttflux(istrip,nlay), stqflux(istrip,nlay)        _RL sttflux(istrip,nlay), stqflux(istrip,nlay)
213        real frqtrb(istrip,nlay-1)        _RL frqtrb(istrip,nlay-1)
214        real dshdthg(istrip,nlay),dthdthg(istrip,nlay)        _RL dshdthg(istrip,nlay),dthdthg(istrip,nlay)
215        real dshdshg(istrip,nlay),dthdshg(istrip,nlay)        _RL dshdshg(istrip,nlay),dthdshg(istrip,nlay)
216        real transth(istrip,nlay), transsh(istrip,nlay)        _RL transth(istrip,nlay), transsh(istrip,nlay)
217        real checktrb(istrip,nlay)  
218          _RL tc(istrip),td(istrip),qa(istrip)
219        real tc(istrip),td(istrip),qa(istrip)        _RL swet1(istrip),swet2(istrip),swet3(istrip)
220        real swet1(istrip),swet2(istrip),swet3(istrip)        _RL capacity(istrip),snowdepth(istrip)
221        real capacity(istrip),snowdepth(istrip)        _RL stz0(istrip)
222        real strflx(istrip), stz0(istrip)        _RL stdiag(istrip)
223        real stcndi(istrip), stdiag(istrip)        _RL tends(istrip),sustar(istrip), sz0(istrip),pbldpth(istrip)
224        real tends(istrip),sustar(istrip), sz0(istrip),pbldpth(istrip)        _RL sct(istrip), scu(istrip), swinds(istrip)
225        real sdtsrf(istrip), stg(istrip), sqs(istrip)        _RL stu2m(istrip),stv2m(istrip),stt2m(istrip),stq2m(istrip)
226        real sct(istrip), scu(istrip), swinds(istrip), sdtg(istrip)        _RL stu10m(istrip),stv10m(istrip),stt10m(istrip),stq10m(istrip)
       real sts(istrip), sqg(istrip)  
       real stu2m(istrip),stv2m(istrip),stt2m(istrip),stq2m(istrip)  
       real stu10m(istrip),stv10m(istrip),stt10m(istrip),stq10m(istrip)  
227        integer  stwatr(istrip)        integer  stwatr(istrip)
228        real  wspeed(istrip)        _RL  wspeed(istrip)
229    
230        real ctsave(istrip),xxsave(istrip),yysave(istrip),zetasave(istrip)        _RL ctsave(istrip),xxsave(istrip),yysave(istrip)
231        real xlsave(istrip,nlay),khsave(istrip,nlay)        _RL zetasave(istrip)
232        real qliq(istrip,nlay),turbfcc(istrip,nlay)        _RL xlsave(istrip,nlay),khsave(istrip,nlay)
233        real qliqmsc(nchp,nlay),fccmsc(nchp,nlay)        _RL qliq(istrip,nlay),turbfcc(istrip,nlay)
234          _RL qliqmsc(nchp,nlay),fccmsc(nchp,nlay)
235    
236        integer ndlsm        integer ndlsm
237        parameter ( ndlsm = 1 )        parameter ( ndlsm = 1)
238        real qdiaglsm(nchp,ndlsm)        _RL qdiaglsm(nchp,ndlsm)
   
       integer nsecf,nmonf,ndayf  
       nsecf(n) = n/10000*3600 + mod(n,10000)/100* 60 + mod(n,100)  
       nmonf(n) = mod(n,10000)/100  
       ndayf(n) = mod(n,100)  
239    
240        real pi,secday,sdayopi2,rgas,akap,cp,alhl        _RL pi,secday,sdayopi2,rgas,akap,cp,alhl
241        real faceps,grav,caltoj,virtcon,getcon        _RL faceps,grav,caltoj,virtcon,getcon
242        real heatw,undef,timstp,delttrb,dttrb,ra        _RL heatw,undef,timstp,delttrb,dttrb,ra
243        real edle,rmu,cltj10,atimstp,tice,const        _RL edle,rmu,cltj10,atimstp,tice,const
244        integer istnp1,istnlay,itrtrb,i,j,L,k,nn,nt        integer istnp1,istnlay,itrtrb,i,j,L,nn,nt
245        integer nocean, nice        integer nocean, nice
246        integer ndmoist,time_left,ndum        integer ndmoist,time_left,ndum
247        integer ntracedim        integer ntracedim
248        real    dtfac,timstp2,sum        _RL    dtfac,timstp2,sum
249  C logical begin flag - set to true to indicate a cold start  C logical begin flag - set to true to indicate a cold start
250        logical qbeg            logical qbeg    
251    
252  #if CRAY        integer n,nsecf,nmonf,ndayf
253  #if f77        nsecf(n) = n/10000*3600 + mod(n,10000)/100* 60 + mod(n,100)
254          nmonf(n) = mod(n,10000)/100
255          ndayf(n) = mod(n,100)
256    
257    #ifdef CRAY
258    #ifdef f77
259  cfpp$ expand (qsat)  cfpp$ expand (qsat)
260  #endif  #endif
261  #endif  #endif
262    
263  c   compute variables that do not change  c   compute variables that do not change
264  c  c
265    
266         pi = 4.*atan(1.)         pi = 4.*atan(1.)
267         secday   = getcon('SDAY')         secday   = getcon('SDAY')
268         sdayopi2 = getcon('SDAY') / (pi*2.)         sdayopi2 = getcon('SDAY') / (pi*2.)
# Line 341  c ************************************** Line 309  c **************************************
309    
310        sum = 0.0        sum = 0.0
311        do L=1,nlay        do L=1,nlay
312        do n=1,nchp        do n=1,nchptot
313        sum = sum + tke(n,L)        sum = sum + tke(n,L)
314        enddo        enddo
315        enddo        enddo
316    
317    #ifdef ALLOW_USE_MPI
318        call mpi_allreduce(sum,const,1,mpi_double_precision,mpi_sum,        call mpi_allreduce(sum,const,1,mpi_double_precision,mpi_sum,
319       .                                                           comm,n)       .                                                mpi_comm_world,n)
320    #else
321          const = sum
322    #endif
323    
324        if( const.eq.0.0 ) then        if( const.eq.0.0 ) then
325        qbeg = .true.        qbeg = .true.
# Line 359  c ************************************** Line 332  c **************************************
332            endif            endif
333        endif        endif
334    
335          if(itground.gt.0) then
336          do j =1,jm
337          do i =1,im
338          qdiag(i,j,itground,bi,bj) = qdiag(i,j,itground,bi,bj) + tgz(i,j)
339          enddo
340          enddo
341          endif
342    
343  c **********************************************************************  c **********************************************************************
344  c                            Initialization  c                            Initialization
345  c **********************************************************************  c **********************************************************************
346                
 c Zero-out 2m and 10m Couplings  
 c -----------------------------  
       do j = 1,jm  
       do i = 1,im  
        u2m(i,j) = 0.0  
        v2m(i,j) = 0.0  
        t2m(i,j) = 0.0  
        q2m(i,j) = 0.0  
       u10m(i,j) = 0.0  
       v10m(i,j) = 0.0  
       t10m(i,j) = 0.0  
       q10m(i,j) = 0.0  
       enddo  
       enddo  
   
347  c Initialize diagnostic for ground temperature change  c Initialize diagnostic for ground temperature change
348  c ---------------------------------------------------  c ---------------------------------------------------
349        if(idtg.gt.0) then        if(idtg.gt.0) then
350        do i =1,ijall        do j =1,jm
351        qdiag(i,1,idtg) = qdiag(i,1,idtg) - tgz(i,1)        do i =1,im
352          qdiag(i,j,idtg,bi,bj) = qdiag(i,j,idtg,bi,bj) - tgz(i,j)
353          enddo
354        enddo        enddo
355        endif        endif
356    
# Line 392  c       do conversion of model state var Line 360  c       do conversion of model state var
360  c        (ocean points appended to tile space land point arrays)  c        (ocean points appended to tile space land point arrays)
361  c **********************************************************************  c **********************************************************************
362    
363        numstrips   = ((nchp-1)/istrip) + 1        numstrips   = ((nchptot-1)/istrip) + 1
364    
365          call grd2msc(pz(1,1),im,jm,igrd,pmsc,nchp,nchptot)
366    
367        call grd2msc(pz(1,1),im,jm,igrd,pmsc,nchp,nchp)        call grd2msc(tgz,im,jm,igrd,tground,nchp,nchptot)
       call grd2msc(tgz,im,jm,igrd,tground,nchp,nchp)  
368        do i = 1,ijall        do i = 1,ijall
369         tmpdiag(i,1) = st4(i,1) + dst4(i,1)*(tgz(i,1)-tempref(i,1))         tmpdiag(i,1) = st4(i,1) + dst4(i,1)*(tgz(i,1)-tempref(i,1))
370       1                         - dst4(i,1)* tgz(i,1)       1                         - dst4(i,1)* tgz(i,1)
371        enddo        enddo
372        call grd2msc(tmpdiag,im,jm,igrd,alwcoeff,nchp,nchp)        call grd2msc(tmpdiag,im,jm,igrd,alwcoeff,nchp,nchptot)
373        do i = 1,ijall        do i = 1,ijall
374         tmpdiag(i,1) = dst4(i,1)         tmpdiag(i,1) = dst4(i,1)
375        enddo        enddo
376        call grd2msc(tmpdiag,im,jm,igrd,blwcoeff,nchp,nchp)        call grd2msc(tmpdiag,im,jm,igrd,blwcoeff,nchp,nchptot)
377        do i = 1,ijall        do i = 1,ijall
378         tmpdiag(i,1) = fdifpar(i,1) * radswt(i,1)         tmpdiag(i,1) = fdifpar(i,1) * radswt(i,1)
379        enddo        enddo
380        call grd2msc(tmpdiag,im,jm,igrd,pardiff,nchp,nchp)        call grd2msc(tmpdiag,im,jm,igrd,pardiff,nchp,nchptot)
381        do i = 1,ijall        do i = 1,ijall
382         tmpdiag(i,1) = fdirpar(i,1) * radswt(i,1)         tmpdiag(i,1) = fdirpar(i,1) * radswt(i,1)
383        enddo        enddo
384        call grd2msc(tmpdiag,im,jm,igrd,pardirct,nchp,nchp)        call grd2msc(tmpdiag,im,jm,igrd,pardirct,nchp,nchptot)
385        do i = 1,ijall        do i = 1,ijall
386         tmpdiag(i,1) = radswg(i,1) * radswt(i,1)         tmpdiag(i,1) = radswg(i,1) * radswt(i,1)
387        enddo        enddo
388        call grd2msc(tmpdiag,im,jm,igrd,netsw,nchp,nchp)        call grd2msc(tmpdiag,im,jm,igrd,netsw,nchp,nchptot)
389        do i = 1,ijall        do i = 1,ijall
390         tmpdiag(i,1) = radlwg(i,1) + dst4(i,1)*(tgz(i,1)-tempref(i,1))         tmpdiag(i,1) = radlwg(i,1) + dst4(i,1)*(tgz(i,1)-tempref(i,1))
391        enddo        enddo
392        call grd2msc(tmpdiag,im,jm,igrd,netlw,nchp,nchp)        call grd2msc(tmpdiag,im,jm,igrd,netlw,nchp,nchptot)
393        call grd2msc(thkz,im,jm,igrd,icethk,nchp,nchp)        call grd2msc(thkz,im,jm,igrd,icethk,nchp,nchptot)
394        call grd2msc(rainlsp,im,jm,igrd,lsprec,nchp,nchp)        call grd2msc(rainlsp,im,jm,igrd,lsprec,nchp,nchptot)
395        call grd2msc(rainconv,im,jm,igrd,cnvprec,nchp,nchp)        call grd2msc(rainconv,im,jm,igrd,cnvprec,nchp,nchptot)
396        call grd2msc(snowfall,im,jm,igrd,snowprec,nchp,nchp)        call grd2msc(snowfall,im,jm,igrd,snowprec,nchp,nchptot)
397    
398  C Call chpprm to get non-varying vegetation and soil characteristics  C Call chpprm to get non-varying vegetation and soil characteristics
399    
# Line 437  c ************************************** Line 406  c **************************************
406    
407  c   set water  c   set water
408    
409        do i = 1,nchp        do i = 1,nchptot
410         water(i) = 0         water(i) = 0
411         if((ityp(i).eq.100).and.(icethk(i).eq.0. ))water(i) = 1         if((ityp(i).eq.100).and.(icethk(i).eq.0. ))water(i) = 1
412        enddo        enddo
413    
414  c   roughness length z0  c   roughness length z0
415  c  c
416        do i =1,nchp        do i =1,nchptot
417         if (icethk(i).gt.0.) then         if (icethk(i).gt.0.) then
418          z0(i) = 1.e-4          z0(i) = 1.e-4
419         else if (ityp(i).eq.100) then         else if (ityp(i).eq.100) then
# Line 463  c  (it has sst from the tgz array over t Line 432  c  (it has sst from the tgz array over t
432    
433  C value of sh at ground  C value of sh at ground
434  C ---------------------  C ---------------------
435        do I =1,nchp        do I =1,nchptot
436        utility(I) = pmsc(i) + ptop        utility(I) = pmsc(i) + ptop
437        call qsat ( tground(i),utility(i),shg(i),dqsdt(i),.true. )        call qsat ( tground(i),utility(i),shg(i),dqsdt(i),.true. )
438        enddo        enddo
# Line 474  c  (it has qstar at tground over the sea Line 443  c  (it has qstar at tground over the sea
443        do i = 1,nchplnd        do i = 1,nchplnd
444         qground(i) = ecanopy(i)         qground(i) = ecanopy(i)
445        enddo        enddo
446        do i = nchplnd+1,nchp        do i = nchplnd+1,nchptot
447         qground(i) = shg(i)         qground(i) = shg(i)
448        enddo        enddo
449    
450  c Fill Array Swetshal with Value 1. over oceans and sea ice  c Fill Array Swetshal with Value 1. over oceans and sea ice
451        do i = nchplnd+1,nchp        do i = nchplnd+1,nchptot
452         swetshal(i) = 1.         swetshal(i) = 1.
453        enddo        enddo
454    
455  c compute heat conduction through ice  c compute heat conduction through ice
456  c -----------------------------------  c -----------------------------------
457        const = ( cti / hice ) * cltj10        const = ( cti / hice ) * cltj10
458        do i =1,nchp        do i =1,nchptot
459               qice(i) =  0.0               qice(i) =  0.0
460              dqice(i) =  0.0              dqice(i) =  0.0
461         if( icethk(i).gt.0.0 ) then         if( icethk(i).gt.0.0 ) then
# Line 496  c ----------------------------------- Line 465  c -----------------------------------
465        enddo        enddo
466    
467        if( iqice.gt.0 ) then        if( iqice.gt.0 ) then
468        tmpdiag(:,:) = 0.0        do i =1,ijall
469        call msc2grd (igrd,chfr,qice,nchp,nchp,fracland,tmpdiag,im,jm)        tmpdiag(i,1) = 0.0
470        qdiag(:,:,iqice) = qdiag(:,:,iqice) + tmpdiag(:,:)        enddo
471          call msc2grd (igrd,chfr,qice,nchp,nchptot,fracland,tmpdiag,im,jm)
472          do j =1,jm
473          do i =1,im
474          qdiag(i,j,iqice,bi,bj) = qdiag(i,j,iqice,bi,bj) + tmpdiag(i,j)
475          enddo
476          enddo
477        nqice = nqice + 1        nqice = nqice + 1
478        endif        endif
479    
# Line 517  c*************************************** Line 492  c***************************************
492         call strip2tile(plze,igrd,pe,nchp,ijall,istrip,nlay+1,nn)         call strip2tile(plze,igrd,pe,nchp,ijall,istrip,nlay+1,nn)
493         call strip2tile(pkz,igrd,pk,nchp,ijall,istrip,nlay,nn)         call strip2tile(pkz,igrd,pk,nchp,ijall,istrip,nlay,nn)
494         call strip2tile(pkht,igrd,pke,nchp,ijall,istrip,nlay+1,nn)         call strip2tile(pkht,igrd,pke,nchp,ijall,istrip,nlay+1,nn)
495         do nt = 1,ntracers-ptracers  c      do nt = 1,ntracers-ptracers
496         call strip2tile(qz(1,1,1,ptracers+nt),igrd,tracers(1,1,nt),nchp,  c      call strip2tile(qz(1,1,1,ptracers+nt),igrd,tracers(1,1,nt),nchp,
497       1                                             ijall,istrip,nlay,nn)  c    1                                             ijall,istrip,nlay,nn)
498         enddo  c      enddo
499    
500         call stripit  (z0,stz0,nchp,nchp,istrip,1,nn)         call stripit  (z0,stz0,nchptot,nchp,istrip,1,nn)
501         call stripit  (tground,th(1,nlay+1),nchp,nchp,istrip,1,nn)         call stripit  (tground,th(1,nlay+1),nchptot,nchp,istrip,1,nn)
502         call stripit  (pmsc,pe(1,nlay+1),nchp,nchp,istrip,1,nn)         call stripit  (pmsc,pe(1,nlay+1),nchptot,nchp,istrip,1,nn)
503         call stripit  (tke,qq,nchp,nchp,istrip,nlay-1,nn)         call stripit  (tke,qq,nchptot,nchp,istrip,nlay-1,nn)
504         call stripit  (ctmt,ctsave,nchp,nchp,istrip,1,nn)         call stripit  (ctmt,ctsave,nchptot,nchp,istrip,1,nn)
505         call stripit  (xxmt,xxsave,nchp,nchp,istrip,1,nn)         call stripit  (xxmt,xxsave,nchptot,nchp,istrip,1,nn)
506         call stripit  (yymt,yysave,nchp,nchp,istrip,1,nn)         call stripit  (yymt,yysave,nchptot,nchp,istrip,1,nn)
507         call stripit  (zetamt,zetasave,nchp,nchp,istrip,1,nn)         call stripit  (zetamt,zetasave,nchptot,nchp,istrip,1,nn)
508         call stripit  (xlmt,xlsave,nchp,nchp,istrip,nlay,nn)         call stripit  (xlmt,xlsave,nchptot,nchp,istrip,nlay,nn)
509         call stripit  (khmt,khsave,nchp,nchp,istrip,nlay,nn)         call stripit  (khmt,khsave,nchptot,nchp,istrip,nlay,nn)
510         call stripitint (water,stwatr,nchp,nchp,istrip,1,nn)         call stripitint (water,stwatr,nchptot,nchp,istrip,1,nn)
511    
512         call stripitint (igrd,igrdstr,nchp,nchp,istrip,1,nn)         call stripitint (igrd,igrdstr,nchptot,nchp,istrip,1,nn)
513         call stripit  (chfr,chfrstr,nchp,nchp,istrip,1,nn)         call stripit  (chfr,chfrstr,nchptot,nchp,istrip,1,nn)
514         call stripit  (icethk,icest,nchp,nchp,istrip,1,nn)         call stripit  (icethk,icest,nchptot,nchp,istrip,1,nn)
515         call stripit  (pardiff,pardf,nchp,nchp,istrip,1,nn)         call stripit  (pardiff,pardf,nchptot,nchp,istrip,1,nn)
516         call stripit  (pardirct,pardr,nchp,nchp,istrip,1,nn)         call stripit  (pardirct,pardr,nchptot,nchp,istrip,1,nn)
517         call stripit  (chlt,lats,nchp,nchp,istrip,1,nn)         call stripit  (chlt,lats,nchptot,nchp,istrip,1,nn)
518         call stripit  (chlon,lons,nchp,nchp,istrip,1,nn)         call stripit  (chlon,lons,nchptot,nchp,istrip,1,nn)
519         call stripit  (lsprec,rainls,nchp,nchp,istrip,1,nn)         call stripit  (lsprec,rainls,nchptot,nchp,istrip,1,nn)
520         call stripit  (cnvprec,raincon,nchp,nchp,istrip,1,nn)         call stripit  (cnvprec,raincon,nchptot,nchp,istrip,1,nn)
521         call stripit  (snowprec,newsnow,nchp,nchp,istrip,1,nn)         call stripit  (snowprec,newsnow,nchptot,nchp,istrip,1,nn)
522         call stripit  (netsw,swnet,nchp,nchp,istrip,1,nn)         call stripit  (netsw,swnet,nchptot,nchp,istrip,1,nn)
523         call stripit  (netlw,lwstrip,nchp,nchp,istrip,1,nn)         call stripit  (netlw,lwstrip,nchptot,nchp,istrip,1,nn)
524         call stripit  (alwcoeff,alwrad,nchp,nchp,istrip,1,nn)         call stripit  (alwcoeff,alwrad,nchptot,nchp,istrip,1,nn)
525         call stripit  (blwcoeff,blwrad,nchp,nchp,istrip,1,nn)         call stripit  (blwcoeff,blwrad,nchptot,nchp,istrip,1,nn)
526         call stripit  (alai,laistrip,nchp,nchp,istrip,1,nn)         call stripit  (alai,laistrip,nchptot,nchp,istrip,1,nn)
527         call stripit  (agrn,grnstrip,nchp,nchp,istrip,1,nn)         call stripit  (agrn,grnstrip,nchptot,nchp,istrip,1,nn)
528         call stripit  (z2ch,z2str,nchp,nchp,istrip,1,nn)         call stripit  (z2ch,z2str,nchptot,nchp,istrip,1,nn)
529         call stripit  (sqscat,scatstr,nchp,nchp,istrip,1,nn)         call stripit  (sqscat,scatstr,nchptot,nchp,istrip,1,nn)
530         call stripit  (rsoil1,rs1str,nchp,nchp,istrip,1,nn)         call stripit  (rsoil1,rs1str,nchptot,nchp,istrip,1,nn)
531         call stripit  (rsoil2,rs2str,nchp,nchp,istrip,1,nn)         call stripit  (rsoil2,rs2str,nchptot,nchp,istrip,1,nn)
532         call stripit  (rdc,rdcstr,nchp,nchp,istrip,1,nn)         call stripit  (rdc,rdcstr,nchptot,nchp,istrip,1,nn)
533         call stripit  (u2fac,u2fstr,nchp,nchp,istrip,1,nn)         call stripit  (u2fac,u2fstr,nchptot,nchp,istrip,1,nn)
534         call stripit  (shg,shgstr,nchp,nchp,istrip,1,nn)         call stripit  (shg,shgstr,nchptot,nchp,istrip,1,nn)
535         call stripit  (dqsdt,dqsdtstr,nchp,nchp,istrip,1,nn)         call stripit  (dqsdt,dqsdtstr,nchptot,nchp,istrip,1,nn)
536         call stripit  ( qice, qicestr,nchp,nchp,istrip,1,nn)         call stripit  ( qice, qicestr,nchptot,nchp,istrip,1,nn)
537         call stripit  (dqice,dqicestr,nchp,nchp,istrip,1,nn)         call stripit  (dqice,dqicestr,nchptot,nchp,istrip,1,nn)
538         call stripitint (ityp,types,nchp,nchp,istrip,1,nn)         call stripitint (ityp,types,nchptot,nchp,istrip,1,nn)
539    
540         call stripit  (tground,tc,nchp,nchp,istrip,1,nn)         call stripit  (tground,tc,nchptot,nchp,istrip,1,nn)
541         call stripit  (tdeep,td,nchp,nchp,istrip,1,nn)         call stripit  (tdeep,td,nchptot,nchp,istrip,1,nn)
542         call stripit  (qground,qa,nchp,nchp,istrip,1,nn)         call stripit  (qground,qa,nchptot,nchp,istrip,1,nn)
543         call stripit  (swetshal,swet1,nchp,nchp,istrip,1,nn)         call stripit  (swetshal,swet1,nchptot,nchp,istrip,1,nn)
544         call stripit  (swetroot,swet2,nchp,nchp,istrip,1,nn)         call stripit  (swetroot,swet2,nchptot,nchp,istrip,1,nn)
545         call stripit  (swetdeep,swet3,nchp,nchp,istrip,1,nn)         call stripit  (swetdeep,swet3,nchptot,nchp,istrip,1,nn)
546         call stripit  (snodep,snowdepth,nchp,nchp,istrip,1,nn)         call stripit  (snodep,snowdepth,nchptot,nchp,istrip,1,nn)
547         call stripit  (capac,capacity,nchp,nchp,istrip,1,nn)         call stripit  (capac,capacity,nchptot,nchp,istrip,1,nn)
548    
549         call astro ( nymd,nhms,lats,lons,istrip,cosz,ra )         call astro ( nymd,nhms,lats,lons,istrip,cosz,ra )
550    
# Line 600  c -------------------------------------- Line 575  c --------------------------------------
575        sh(i,nlay+1) = qa(i)        sh(i,nlay+1) = qa(i)
576        enddo        enddo
577    
578        if(iqg.gt.0) then  c     if(iqg.gt.0) then
579        do i=1,istrip  c     do i=1,istrip
580        tmpstrip(i) = sh(i,nlay+1)*1000  c     tmpstrip(i) = sh(i,nlay+1)*1000
581        enddo  c     enddo
582        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
583       1                        qdiag(1,1,iqg),ijall,1,nn,.false.)  c    1                        qdiag(1,1,iqg,bi,bj),ijall,1,nn,.false.)
584        endif  c     endif
585    
586  c value of tracers at the ground  c value of tracers at the ground
587  c ------------------------------  c ------------------------------
588        do nt = 1,ntracers-ptracers  c     do nt = 1,ntracers-ptracers
589         do i = 1,istrip  C      do i = 1,istrip
590          tracers(i,nlay+1,nt) = 0.  C       tracers(i,nlay+1,nt) = 0.
591         enddo  C      enddo
592        enddo  C     enddo
593    
594  c compute virtual potential temperatures  c compute virtual potential temperatures
595  c --------------------------------------  c --------------------------------------
# Line 650  c ------------------------------- Line 625  c -------------------------------
625    
626  c increment diagnostic arrays for quantities calculated before trbfl  c increment diagnostic arrays for quantities calculated before trbfl
627  c ------------------------------------------------------------------  c ------------------------------------------------------------------
628        do i =1,istrip  c     do i =1,istrip
629        stdiag(i) = ( thv(i,nlay+1)-thv(i,nlay) ) / pke(i,nlay+1)  c     stdiag(i) = ( thv(i,nlay+1)-thv(i,nlay) ) / pke(i,nlay+1)
630        enddo  c     enddo
631        if(idtsrf.gt.0) then  c     if(idtsrf.gt.0) then
632        call paste2grd(stdiag,igrd,chfrstr,istrip,nchp,  c     call paste2grd(stdiag,igrd,chfrstr,istrip,nchp,
633       1                      qdiag(1,1,idtsrf),ijall,1,nn,.false.)  c    1                      qdiag(1,1,idtsrf,bi,bj),ijall,1,nn,.false.)
634    c     endif
635    
636    
637          if(1.eq.1)then
638          print *,' In turb before trbflx - strip ',nn,' out of ',numstrips
639          print *,' bi = ',bi
640          print *,' ntracers ',ntracers,' ptracers ',ptracers
641          print *,'dttrb,itrtrb,rmu,edle ',dttrb,' ',itrtrb,' ',rmu,' ',edle
642          print *,' nchp ',nchp,' nchptot ',nchptot,' nchplnd ',nchplnd
643          print *,' qbeg, tprof ',qbeg,' ',tprof
644          print *,'istrip,nlay,nymd,nhms ',istrip,' ',nlay,' ',nymd,' ',nhms
645          print *,' grav,cp,rgas,faceps,virtcon,undef ',
646         .     grav,' ',cp,' ',rgas,' ',faceps,' ',virtcon,' ',undef
647          print *,' field: th ',th
648    c     print *,' field: thv ',thv
649    c     print *,' field: sh ',sh
650    c     print *,' field: u ',u
651    c     print *,' field: v ',v
652          print *,' field: p ',p
653    c     print *,' field: pe ',pe
654    c     print *,' field: pk ',pk
655    c     print *,' field: pke ',pke
656    c     print *,' field: dpstr ',dpstr
657    c     print *,' field: stwatr ',stwatr
658    c     print *,' field: stz0 ',stz0
659          endif
660    
661          if(iudiag1.gt.0) then
662          call paste2grd(th,igrd,chfrstr,istrip,nchp,
663         1                   tempor1,ijall,nlay,nn,.false.)
664          endif
665          if(iudiag2.gt.0) then
666          call paste2grd(p,igrd,chfrstr,istrip,nchp,
667         1                   tempor2,ijall,nlay,nn,.false.)
668        endif        endif
669    
670    
671  c call trbflx  c call trbflx
672  c -----------  c -----------
673        call trbflx(nn,th,thv,sh,u,v,qq,p,pe,pk,pke,dpstr,stwatr,stz0,        call trbflx(nn,th,thv,sh,u,v,qq,p,pe,pk,pke,dpstr,stwatr,stz0,
# Line 667  c ----------- Line 677  c -----------
677       4 stq10m,istrip,nlay,nymd,nhms,grav,cp,rgas,faceps,virtcon,undef,       4 stq10m,istrip,nlay,nymd,nhms,grav,cp,rgas,faceps,virtcon,undef,
678       5 dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,       5 dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,
679       6 hsturb,dhsdqa,dhsdtc,transth,transsh,       6 hsturb,dhsdqa,dhsdtc,transth,transsh,
680       7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc,       7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)
681       8 checktrb)  
682    
683          if(iudiag3.gt.0) call paste2grd(th,igrd,chfrstr,istrip,nchp,
684         1                   tempor3,ijall,nlay,nn,.false.)
685    
686        call pastit (qq,tke,istrip,nchp,nchp,nlay,nn)        if(2.eq.1)then
687        call pastit (ctsave,ctmt,istrip,nchp,nchp,1,nn)        print *,' In turbio, Just after trbflx for strip ',nn,' bi = ',bi
688        call pastit (xxsave,xxmt,istrip,nchp,nchp,1,nn)        print *,' field: stuflux ',stuflux
689        call pastit (yysave,yymt,istrip,nchp,nchp,1,nn)        print *,' field: stvflux ',stvflux
690        call pastit (zetasave,zetamt,istrip,nchp,nchp,1,nn)        print *,' field: dshdthg ',dshdthg
691        call pastit (xlsave,xlmt,istrip,nchp,nchp,nlay,nn)        print *,' field: dshdshg ',dshdshg
692        call pastit (khsave,khmt,istrip,nchp,nchp,nlay,nn)        print *,' field: dthdthg ',dthdthg
693          print *,' field: dthdshg ',dthdshg
694          print *,' field: scu ',scu
695          print *,' field: eturb ',eturb
696          print *,' field: dedqa ',dedqa
697          print *,' field: dedtc ',dedtc
698          print *,' field: hsturb ',hsturb
699          print *,' field: dhsdqa ',dhsdqa
700          print *,' field: dhsdtc ',dhsdtc
701          print *,' field: transth ',transth
702          print *,' field: transsh ',transsh
703          endif
704    
705          call pastit (qq,tke,istrip,nchp,nchptot,nlay,nn)
706          call pastit (ctsave,ctmt,istrip,nchp,nchptot,1,nn)
707          call pastit (xxsave,xxmt,istrip,nchp,nchptot,1,nn)
708          call pastit (yysave,yymt,istrip,nchp,nchptot,1,nn)
709          call pastit (zetasave,zetamt,istrip,nchp,nchptot,1,nn)
710          call pastit (xlsave,xlmt,istrip,nchp,nchptot,nlay,nn)
711          call pastit (khsave,khmt,istrip,nchp,nchptot,nlay,nn)
712    
713        call pastit (qliq  ,qliqmsc,istrip,nchp,nchp,nlay,nn)        call pastit (qliq  ,qliqmsc,istrip,nchp,nchptot,nlay,nn)
714        call pastit (turbfcc,fccmsc,istrip,nchp,nchp,nlay,nn)        call pastit (turbfcc,fccmsc,istrip,nchp,nchptot,nlay,nn)
715    
716  c  New diagnostic: potential evapotranspiration  c  New diagnostic: potential evapotranspiration
717        do i = 1,istrip        do i = 1,istrip
# Line 701  C*************************************** Line 733  C***************************************
733         hlwdwn(i) = alwrad(i)+blwrad(i)*tc(i)-lwstrip(i)         hlwdwn(i) = alwrad(i)+blwrad(i)*tc(i)-lwstrip(i)
734         psurf(i) = pe(i,nlay+1)         psurf(i) = pe(i,nlay+1)
735         wspeed(i) = sqrt(u(i,nlay)*u(i,nlay) + v(i,nlay)*v(i,nlay))         wspeed(i) = sqrt(u(i,nlay)*u(i,nlay) + v(i,nlay)*v(i,nlay))
736           if(wspeed(i) .lt. 1.e-20) wspeed(i) = 1.e-20
737  C Note: This LSM precip bug needs to be cleaned up  C Note: This LSM precip bug needs to be cleaned up
738  ccc   newsnow(i) = newsnow(i)*dtfac    ccc   newsnow(i) = newsnow(i)*dtfac  
739  ccc   raincon(i) = raincon(i)*dtfac    ccc   raincon(i) = raincon(i)*dtfac  
# Line 740  ccc   rainls (i) = rainls (i)*dtfac Line 773  ccc   rainls (i) = rainls (i)*dtfac
773  c**********************************************************************  c**********************************************************************
774  c   diagnostics: fill arrays for lsm input fields  c   diagnostics: fill arrays for lsm input fields
775  c**********************************************************************  c**********************************************************************
776        if(isnowfall.gt.0) then  c     if(isnowfall.gt.0) then
777        do i = 1,istrip  c     do i = 1,istrip
778        tmpstrip(i) = newsnow(i)*86400    c     tmpstrip(i) = newsnow(i)*86400  
779        enddo  c     enddo
780        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
781       1                       qdiag(1,1,isnowfall),ijall,1,nn,.false.)  c    1                   qdiag(1,1,isnowfall,bi,bj),ijall,1,nn,.false.)
782        endif  c     endif
783        if(iraincon.gt.0) then  c     if(iraincon.gt.0) then
784        do i = 1,istrip  c     do i = 1,istrip
785        tmpstrip(i) = raincon(i)*86400    c     tmpstrip(i) = raincon(i)*86400  
786        enddo  c     enddo
787        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
788       1                       qdiag(1,1,iraincon),ijall,1,nn,.false.)  c    1                   qdiag(1,1,iraincon,bi,bj),ijall,1,nn,.false.)
789        endif  c     endif
790        if(irainlsp.gt.0) then  c     if(irainlsp.gt.0) then
791        do i = 1,istrip  c     do i = 1,istrip
792        tmpstrip(i) = rainls(i)*86400    c     tmpstrip(i) = rainls(i)*86400  
793        enddo  c     enddo
794        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
795       1                      qdiag(1,1,irainlsp),ijall,1,nn,.false.)  c    1                  qdiag(1,1,irainlsp,bi,bj),ijall,1,nn,.false.)
796        endif  c     endif
797        if(igreen.gt.0) then  c     if(igreen.gt.0) then
798        call paste2grd(grnstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(grnstrip,igrd,chfrstr,istrip,nchp,
799       1                        qdiag(1,1,igreen),ijall,1,nn,.false.)  c    1                    qdiag(1,1,igreen,bi,bj),ijall,1,nn,.false.)
800        endif  c     endif
801        if(ilai.gt.0) then  c     if(ilai.gt.0) then
802        call paste2grd(laistrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(laistrip,igrd,chfrstr,istrip,nchp,
803       1                        qdiag(1,1,ilai),ijall,1,nn,.false.)  c    1                    qdiag(1,1,ilai,bi,bj),ijall,1,nn,.false.)
804        endif  c     endif
805        if(ipardr.gt.0) then  c     if(ipardr.gt.0) then
806        call paste2grd(pardr,igrd,chfrstr,istrip,nchp,  c     call paste2grd(pardr,igrd,chfrstr,istrip,nchp,
807       1                     qdiag(1,1,ipardr),ijall,1,nn,.false.)  c    1                 qdiag(1,1,ipardr,bi,bj),ijall,1,nn,.false.)
808        endif  c     endif
809        if(ipardf.gt.0) then  c     if(ipardf.gt.0) then
810        call paste2grd(pardf,igrd,chfrstr,istrip,nchp,  c     call paste2grd(pardf,igrd,chfrstr,istrip,nchp,
811       1                     qdiag(1,1,ipardf),ijall,1,nn,.false.)  c    1                 qdiag(1,1,ipardf,bi,bj),ijall,1,nn,.false.)
812        endif  c     endif
813        if(idlwdtc.gt.0) then  c     if(idlwdtc.gt.0) then
814        call paste2grd(blwrad,igrd,chfrstr,istrip,nchp,  c     call paste2grd(blwrad,igrd,chfrstr,istrip,nchp,
815       1                      qdiag(1,1,idlwdtc),ijall,1,nn,.false.)  c    1                  qdiag(1,1,idlwdtc,bi,bj),ijall,1,nn,.false.)
816        endif  c     endif
817        if(idhdtc.gt.0) then  c     if(idhdtc.gt.0) then
818        call paste2grd(dhsdtc,igrd,chfrstr,istrip,nchp,  c     call paste2grd(dhsdtc,igrd,chfrstr,istrip,nchp,
819       1                      qdiag(1,1,idhdtc),ijall,1,nn,.false.)  c    1                  qdiag(1,1,idhdtc,bi,bj),ijall,1,nn,.false.)
820        endif  c     endif
821        if(idedtc.gt.0) then  c     if(idedtc.gt.0) then
822        call paste2grd(dedtc,igrd,chfrstr,istrip,nchp,  c     call paste2grd(dedtc,igrd,chfrstr,istrip,nchp,
823       1                     qdiag(1,1,idedtc),ijall,1,nn,.false.)  c    1                 qdiag(1,1,idedtc,bi,bj),ijall,1,nn,.false.)
824        endif  c     endif
825        if(idhdqa.gt.0) then  c     if(idhdqa.gt.0) then
826        call paste2grd(dhsdqa,igrd,chfrstr,istrip,nchp,  c     call paste2grd(dhsdqa,igrd,chfrstr,istrip,nchp,
827       1                      qdiag(1,1,idhdqa),ijall,1,nn,.false.)  c    1                  qdiag(1,1,idhdqa,bi,bj),ijall,1,nn,.false.)
828        endif  c     endif
829        if(idedqa.gt.0) then  c     if(idedqa.gt.0) then
830        call paste2grd(dedqa,igrd,chfrstr,istrip,nchp,  c     call paste2grd(dedqa,igrd,chfrstr,istrip,nchp,
831       1                     qdiag(1,1,idedqa),ijall,1,nn,.false.)  c    1                 qdiag(1,1,idedqa,bi,bj),ijall,1,nn,.false.)
832        endif  c     endif
833        if(ilwgdown.gt.0) then  c     if(ilwgdown.gt.0) then
834        call paste2grd(hlwdwn,igrd,chfrstr,istrip,nchp,  c     call paste2grd(hlwdwn,igrd,chfrstr,istrip,nchp,
835       1                      qdiag(1,1,ilwgdown),ijall,1,nn,.false.)  c    1                  qdiag(1,1,ilwgdown,bi,bj),ijall,1,nn,.false.)
836        endif  c     endif
837  c**********************************************************************  c**********************************************************************
838    
839        if(nland.gt.0)then        if(nland.gt.0)then
840    
841          if(1.eq.1)then
842          print *,' In turbio, Just before tile for strip ',nn,' bi = ',bi
843          print *,' calling tile for ',nland,' land points '
844          print *,' field: types ',types
845          print *,' field: chfrstr ',chfrstr
846    c     print *,' field: rainls ',rainls
847    c     print *,' field: newsnow ',newsnow
848    c     print *,' field: wspeed ',wspeed
849    c     print *,' field: eturb ',eturb
850    c     print *,' field: dedqa ',dedqa
851    c     print *,' field: dedtc ',dedtc
852    c     print *,' field: hsturb ',hsturb
853    c     print *,' field: dhsdqa ',dhsdqa
854    c     print *,' field: dhsdtc ',dhsdtc
855    c     print *,' field: tmpnlay ',tmpnlay
856    c     print *,' field: sh(nlay) ',(sh(i,nlay),i=1,istrip)
857    c     print *,' field: cd ',cd
858    c     print *,' field: cosz ',cosz
859    c     print *,' field: pardr ',pardr
860    c     print *,' field: pardf ',pardf
861    c     print *,' field: swnet ',swnet
862    c     print *,' field: hlwdwn ',hlwdwn
863    c     print *,' field: psurf ',psurf
864    c     print *,' field: laistrip ',laistrip
865    c     print *,' field: grnstrip ',grnstrip
866    c     print *,' field: z2str ',z2str
867    c     print *,' field: scatstr ',scatstr
868    c     print *,' field: z2str ',z2str
869    c     print *,' field: rs1str ',rs1str
870    c     print *,' field: rs1str ',rs2str
871    c     print *,' field: rdcstr ',rdcstr
872    c     print *,' field: u2fstr ',u2fstr
873    c     print *,' field: dqsdtstr ',dqsdtstr
874    c     print *,' field: alwrad ',alwrad
875    c     print *,' field: blwrad ',blwrad
876          print *,' field: tc ',tc
877    c     print *,' field: td ',td
878    c     print *,' field: swet1 ',swet1
879    c     print *,' field: swet2 ',swet2
880    c     print *,' field: swet3 ',swet3
881    c     print *,' field: capacity ',capacity
882    c     print *,' field: snowdepth ',snowdepth
883          endif
884    
885          if(isdiag1.gt.0) then
886          call paste2grd(tc,igrd,chfrstr,istrip,nchp,
887         1                    tempsfc1,ijall,1,nn,.false.)
888          endif
889    
890    
891         call tile (         call tile (
892       I   nland, timstp, types, rainls, raincon,  newsnow,  wspeed,       I   nland, timstp, types, rainls, raincon,  newsnow,  wspeed,
893       I   eturb,  dedqa,  dedtc,  hsturb, dhsdqa, dhsdtc,       I   eturb,  dedqa,  dedtc,  hsturb, dhsdqa, dhsdtc,
# Line 819  c*************************************** Line 903  c***************************************
903       O   strdg5, strdg6, strdg7, strdg8, strdg9)       O   strdg5, strdg6, strdg7, strdg8, strdg9)
904        endif        endif
905    
906    
907          if(isdiag2.gt.0) then
908          call paste2grd(tc,igrd,chfrstr,istrip,nchp,
909         1                    tempsfc2,ijall,1,nn,.false.)
910          endif
911    
912          if(2.eq.1)then
913          print *,' In turbio, Just after tile for strip ',nn
914          print *,' calling tile for ',nland,' land points '
915          print *,' field: tc ',tc
916          print *,' field: td ',td
917          print *,' field: strdg1 ',strdg1
918          print *,' field: strdg2 ',strdg2
919          print *,' field: strdg3 ',strdg3
920          print *,' field: strdg4 ',strdg4
921          print *,' field: strdg5 ',strdg5
922          print *,' field: strdg6 ',strdg6
923          print *,' field: strdg7 ',strdg7
924          print *,' field: strdg8 ',strdg8
925          print *,' field: strdg9 ',strdg9
926          print *,' field: swet1 ',swet1
927          print *,' field: swet2 ',swet2
928          print *,' field: swet3 ',swet3
929          print *,' field: capacity ',capacity
930          print *,' field: snowdepth ',snowdepth
931          endif
932        if( nice.gt.0 ) then        if( nice.gt.0 ) then
933           print *,' Calling sea ice routine - SHOULD NOT BE HERE!'
934         call seaice ( nocean, timstp,     hice,         call seaice ( nocean, timstp,     hice,
935       .               eturb(nland+1),    dedtc(nland+1),       .               eturb(nland+1),    dedtc(nland+1),
936       .              hsturb(nland+1),   dhsdtc(nland+1),       .              hsturb(nland+1),   dhsdtc(nland+1),
# Line 833  c*************************************** Line 944  c***************************************
944  c   diagnostics: fill arrays for lsm output fields  c   diagnostics: fill arrays for lsm output fields
945  c***********************************************************************  c***********************************************************************
946    
947        if(irunoff.gt.0) then  c     if(irunoff.gt.0) then
948        call paste2grd(runoff,igrd,chfrstr,istrip,nchp,  c     call paste2grd(runoff,igrd,chfrstr,istrip,nchp,
949       1                      qdiag(1,1,irunoff),ijall,1,nn,.false.)  c    1                  qdiag(1,1,irunoff,bi,bj),ijall,1,nn,.false.)
950        endif  c     endif
951        if(ifwsoil.gt.0) then  c     if(ifwsoil.gt.0) then
952        call paste2grd(fwsoil,igrd,chfrstr,istrip,nchp,  c     call paste2grd(fwsoil,igrd,chfrstr,istrip,nchp,
953       1                      qdiag(1,1,ifwsoil),ijall,1,nn,.false.)  c    1                  qdiag(1,1,ifwsoil,bi,bj),ijall,1,nn,.false.)
954        endif  c     endif
955        if(igdrain.gt.0) then  c     if(igdrain.gt.0) then
956        call paste2grd(gdrain,igrd,chfrstr,istrip,nchp,  c     call paste2grd(gdrain,igrd,chfrstr,istrip,nchp,
957       1                      qdiag(1,1,igdrain),ijall,1,nn,.false.)  c    1                  qdiag(1,1,igdrain,bi,bj),ijall,1,nn,.false.)
958        endif  c     endif
959        if(isnowmelt.gt.0) then  c     if(isnowmelt.gt.0) then
960        call paste2grd(smelt,igrd,chfrstr,istrip,nchp,  c     call paste2grd(smelt,igrd,chfrstr,istrip,nchp,
961       1                     qdiag(1,1,isnowmelt),ijall,1,nn,.false.)  c    1                 qdiag(1,1,isnowmelt,bi,bj),ijall,1,nn,.false.)
962        endif  c     endif
963        if(ieveg.gt.0) then  c     if(ieveg.gt.0) then
964        call paste2grd(eveg,igrd,chfrstr,istrip,nchp,  c     call paste2grd(eveg,igrd,chfrstr,istrip,nchp,
965       1                    qdiag(1,1,ieveg),ijall,1,nn,.false.)  c    1                    qdiag(1,1,ieveg,bi,bj),ijall,1,nn,.false.)
966        endif  c     endif
967        if(iesnow.gt.0) then  c     if(iesnow.gt.0) then
968        call paste2grd(esno,igrd,chfrstr,istrip,nchp,  c     call paste2grd(esno,igrd,chfrstr,istrip,nchp,
969       1                    qdiag(1,1,iesnow),ijall,1,nn,.false.)  c    1                    qdiag(1,1,iesnow,bi,bj),ijall,1,nn,.false.)
970        endif  c     endif
971        if(iesoil.gt.0) then  c     if(iesoil.gt.0) then
972        call paste2grd(esoi,igrd,chfrstr,istrip,nchp,  c     call paste2grd(esoi,igrd,chfrstr,istrip,nchp,
973       1                    qdiag(1,1,iesoil),ijall,1,nn,.false.)  c    1                    qdiag(1,1,iesoil,bi,bj),ijall,1,nn,.false.)
974        endif  c     endif
975        if(ieresv.gt.0) then  c     if(ieresv.gt.0) then
976        call paste2grd(eint,igrd,chfrstr,istrip,nchp,  c     call paste2grd(eint,igrd,chfrstr,istrip,nchp,
977       1                    qdiag(1,1,ieresv),ijall,1,nn,.false.)  c    1                    qdiag(1,1,ieresv,bi,bj),ijall,1,nn,.false.)
978        endif  c     endif
979        if(ievpot.gt.0) then  c     if(ievpot.gt.0) then
980        call paste2grd(evpot,igrd,chfrstr,istrip,nchp,  c     call paste2grd(evpot,igrd,chfrstr,istrip,nchp,
981       1                     qdiag(1,1,ievpot),ijall,1,nn,.false.)  c    1                     qdiag(1,1,ievpot,bi,bj),ijall,1,nn,.false.)
982        endif  c     endif
983        if(idtc.gt.0) then  c     if(idtc.gt.0) then
984        call paste2grd(strdg1,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg1,igrd,chfrstr,istrip,nchp,
985       1                      qdiag(1,1,idtc),ijall,1,nn,.false.)  c    1                      qdiag(1,1,idtc,bi,bj),ijall,1,nn,.false.)
986        endif  c     endif
987        if(idqc.gt.0) then  c     if(idqc.gt.0) then
988        call paste2grd(strdg2,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg2,igrd,chfrstr,istrip,nchp,
989       1                      qdiag(1,1,idqc),ijall,1,nn,.false.)  c    1                      qdiag(1,1,idqc,bi,bj),ijall,1,nn,.false.)
990        endif  c     endif
991        if(itcdtc.gt.0) then  c     if(itcdtc.gt.0) then
992        call paste2grd(strdg3,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg3,igrd,chfrstr,istrip,nchp,
993       1                      qdiag(1,1,itcdtc),ijall,1,nn,.false.)  c    1                      qdiag(1,1,itcdtc,bi,bj),ijall,1,nn,.false.)
994        endif  c     endif
995        if(iraddtc.gt.0) then  c     if(iraddtc.gt.0) then
996        call paste2grd(strdg4,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg4,igrd,chfrstr,istrip,nchp,
997       1                      qdiag(1,1,iraddtc),ijall,1,nn,.false.)  c    1                      qdiag(1,1,iraddtc,bi,bj),ijall,1,nn,.false.)
998        endif  c     endif
999        if(isensdtc.gt.0) then  c     if(isensdtc.gt.0) then
1000        call paste2grd(strdg5,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg5,igrd,chfrstr,istrip,nchp,
1001       1                      qdiag(1,1,isensdtc),ijall,1,nn,.false.)  c    1                     qdiag(1,1,isensdtc,bi,bj),ijall,1,nn,.false.)
1002        endif  c     endif
1003        if(ilatdtc.gt.0) then  c     if(ilatdtc.gt.0) then
1004        call paste2grd(strdg6,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg6,igrd,chfrstr,istrip,nchp,
1005       1                      qdiag(1,1,ilatdtc),ijall,1,nn,.false.)  c    1                      qdiag(1,1,ilatdtc,bi,bj),ijall,1,nn,.false.)
1006        endif  c     endif
1007        if(itddtc.gt.0) then  c     if(itddtc.gt.0) then
1008        call paste2grd(strdg7,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg7,igrd,chfrstr,istrip,nchp,
1009       1                      qdiag(1,1,itddtc),ijall,1,nn,.false.)  c    1                      qdiag(1,1,itddtc,bi,bj),ijall,1,nn,.false.)
1010        endif  c     endif
1011        if(iqcdtc.gt.0) then  c     if(iqcdtc.gt.0) then
1012        call paste2grd(strdg8,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg8,igrd,chfrstr,istrip,nchp,
1013       1                      qdiag(1,1,iqcdtc),ijall,1,nn,.false.)  c    1                      qdiag(1,1,iqcdtc,bi,bj),ijall,1,nn,.false.)
1014        endif  c     endif
1015  c**********************************************************************  c***********************************************************************
1016    
1017        if( ndlsm.gt.1 ) then        if( ndlsm.gt.1 ) then
1018        call pstbitint(types,qdiaglsm(1,1),istrip,nchp,nchp,1,nn)        call pstbitint(types,qdiaglsm(1,1),istrip,nchp,nchptot,1,nn)
1019        call pstbmpit(chfrstr,qdiaglsm(1,2),istrip,nchp,nchp,1,nn)        call pstbmpit(chfrstr,qdiaglsm(1,2),istrip,nchp,nchptot,1,nn)
1020        call pstbmpit(lats,qdiaglsm(1,3),istrip,nchp,nchp,1,nn)        call pstbmpit(lats,qdiaglsm(1,3),istrip,nchp,nchptot,1,nn)
1021        call pstbmpit(lons,qdiaglsm(1,4),istrip,nchp,nchp,1,nn)        call pstbmpit(lons,qdiaglsm(1,4),istrip,nchp,nchptot,1,nn)
1022        call pstbmpit(igrdstr,qdiaglsm(1,5),istrip,nchp,nchp,1,nn)  c     call pstbmpit(igrdstr,qdiaglsm(1,5),istrip,nchp,nchptot,1,nn)
1023        call pstbmpit(tc,qdiaglsm(1,6),istrip,nchp,nchp,1,nn)        call pstbmpit(tc,qdiaglsm(1,6),istrip,nchp,nchptot,1,nn)
1024        call pstbmpit(td,qdiaglsm(1,7),istrip,nchp,nchp,1,nn)        call pstbmpit(td,qdiaglsm(1,7),istrip,nchp,nchptot,1,nn)
1025        call pstbmpit(qa,qdiaglsm(1,8),istrip,nchp,nchp,1,nn)        call pstbmpit(qa,qdiaglsm(1,8),istrip,nchp,nchptot,1,nn)
1026        call pstbmpit(swet1,qdiaglsm(1,9),istrip,nchp,nchp,1,nn)        call pstbmpit(swet1,qdiaglsm(1,9),istrip,nchp,nchptot,1,nn)
1027        call pstbmpit(swet2,qdiaglsm(1,10),istrip,nchp,nchp,1,nn)        call pstbmpit(swet2,qdiaglsm(1,10),istrip,nchp,nchptot,1,nn)
1028        call pstbmpit(swet3,qdiaglsm(1,11),istrip,nchp,nchp,1,nn)        call pstbmpit(swet3,qdiaglsm(1,11),istrip,nchp,nchptot,1,nn)
1029        call pstbmpit(capacity,qdiaglsm(1,12),istrip,nchp,nchp,1,nn)        call pstbmpit(capacity,qdiaglsm(1,12),istrip,nchp,nchptot,1,nn)
1030        call pstbmpit(snowdepth,qdiaglsm(1,13),istrip,nchp,nchp,1,nn)        call pstbmpit(snowdepth,qdiaglsm(1,13),istrip,nchp,nchptot,1,nn)
1031        call pstbmpit(eturb,qdiaglsm(1,14),istrip,nchp,nchp,1,nn)        call pstbmpit(eturb,qdiaglsm(1,14),istrip,nchp,nchptot,1,nn)
1032        call pstbmpit(hsturb,qdiaglsm(1,15),istrip,nchp,nchp,1,nn)        call pstbmpit(hsturb,qdiaglsm(1,15),istrip,nchp,nchptot,1,nn)
1033        call pstbmpit(cd,qdiaglsm(1,16),istrip,nchp,nchp,1,nn)        call pstbmpit(cd,qdiaglsm(1,16),istrip,nchp,nchptot,1,nn)
1034        call pstbmpit(laistrip,qdiaglsm(1,17),istrip,nchp,nchp,1,nn)        call pstbmpit(laistrip,qdiaglsm(1,17),istrip,nchp,nchptot,1,nn)
1035        call pstbmpit(grnstrip,qdiaglsm(1,18),istrip,nchp,nchp,1,nn)        call pstbmpit(grnstrip,qdiaglsm(1,18),istrip,nchp,nchptot,1,nn)
1036        call pstbmpit(eint,qdiaglsm(1,19),istrip,nchp,nchp,1,nn)        call pstbmpit(eint,qdiaglsm(1,19),istrip,nchp,nchptot,1,nn)
1037        call pstbmpit(esoi,qdiaglsm(1,20),istrip,nchp,nchp,1,nn)        call pstbmpit(esoi,qdiaglsm(1,20),istrip,nchp,nchptot,1,nn)
1038        call pstbmpit(eveg,qdiaglsm(1,21),istrip,nchp,nchp,1,nn)        call pstbmpit(eveg,qdiaglsm(1,21),istrip,nchp,nchptot,1,nn)
1039        call pstbmpit(esno,qdiaglsm(1,22),istrip,nchp,nchp,1,nn)        call pstbmpit(esno,qdiaglsm(1,22),istrip,nchp,nchptot,1,nn)
1040        call pstbmpit(strdg1,qdiaglsm(1,23),istrip,nchp,nchp,1,nn)        call pstbmpit(strdg1,qdiaglsm(1,23),istrip,nchp,nchptot,1,nn)
1041        call pstbmpit(strdg2,qdiaglsm(1,24),istrip,nchp,nchp,1,nn)        call pstbmpit(strdg2,qdiaglsm(1,24),istrip,nchp,nchptot,1,nn)
1042        call pstbmpit(strdg3,qdiaglsm(1,25),istrip,nchp,nchp,1,nn)        call pstbmpit(strdg3,qdiaglsm(1,25),istrip,nchp,nchptot,1,nn)
1043        call pstbmpit(strdg4,qdiaglsm(1,26),istrip,nchp,nchp,1,nn)        call pstbmpit(strdg4,qdiaglsm(1,26),istrip,nchp,nchptot,1,nn)
1044        call pstbmpit(strdg5,qdiaglsm(1,27),istrip,nchp,nchp,1,nn)        call pstbmpit(strdg5,qdiaglsm(1,27),istrip,nchp,nchptot,1,nn)
1045        call pstbmpit(strdg6,qdiaglsm(1,28),istrip,nchp,nchp,1,nn)        call pstbmpit(strdg6,qdiaglsm(1,28),istrip,nchp,nchptot,1,nn)
1046        call pstbmpit(strdg7,qdiaglsm(1,29),istrip,nchp,nchp,1,nn)        call pstbmpit(strdg7,qdiaglsm(1,29),istrip,nchp,nchptot,1,nn)
1047        call pstbmpit(strdg8,qdiaglsm(1,30),istrip,nchp,nchp,1,nn)        call pstbmpit(strdg8,qdiaglsm(1,30),istrip,nchp,nchptot,1,nn)
1048        call pstbmpit(strdg9,qdiaglsm(1,31),istrip,nchp,nchp,1,nn)        call pstbmpit(strdg9,qdiaglsm(1,31),istrip,nchp,nchptot,1,nn)
1049        call pstbmpit(smelt,qdiaglsm(1,32),istrip,nchp,nchp,1,nn)        call pstbmpit(smelt,qdiaglsm(1,32),istrip,nchp,nchptot,1,nn)
1050        call pstbmpit(gdrain,qdiaglsm(1,33),istrip,nchp,nchp,1,nn)        call pstbmpit(gdrain,qdiaglsm(1,33),istrip,nchp,nchptot,1,nn)
1051        call pstbmpit(runsrf,qdiaglsm(1,34),istrip,nchp,nchp,1,nn)        call pstbmpit(runsrf,qdiaglsm(1,34),istrip,nchp,nchptot,1,nn)
1052        call pstbmpit(fwsoil,qdiaglsm(1,35),istrip,nchp,nchp,1,nn)        call pstbmpit(fwsoil,qdiaglsm(1,35),istrip,nchp,nchptot,1,nn)
1053        call pstbmpit(evpot,qdiaglsm(1,36),istrip,nchp,nchp,1,nn)        call pstbmpit(evpot,qdiaglsm(1,36),istrip,nchp,nchptot,1,nn)
1054        call pstbmpit(stt2m,qdiaglsm(1,37),istrip,nchp,nchp,1,nn)        call pstbmpit(stt2m,qdiaglsm(1,37),istrip,nchp,nchptot,1,nn)
1055        call pstbmpit(stq2m,qdiaglsm(1,38),istrip,nchp,nchp,1,nn)        call pstbmpit(stq2m,qdiaglsm(1,38),istrip,nchp,nchptot,1,nn)
1056        endif        endif
1057    
1058        call pastit (tc,tground,istrip,nchp,nchp,1,nn)        call pastit (tc,tground,istrip,nchp,nchptot,1,nn)
1059        call pastit (td,tdeep,istrip,nchp,nchp,1,nn)        call pastit (td,tdeep,istrip,nchp,nchptot,1,nn)
1060        call pastit (qa,qground,istrip,nchp,nchp,1,nn)        call pastit (qa,qground,istrip,nchp,nchptot,1,nn)
1061        call pastit (swet1,swetshal,istrip,nchp,nchp,1,nn)        call pastit (swet1,swetshal,istrip,nchp,nchptot,1,nn)
1062        call pastit (swet2,swetroot,istrip,nchp,nchp,1,nn)        call pastit (swet2,swetroot,istrip,nchp,nchptot,1,nn)
1063        call pastit (swet3,swetdeep,istrip,nchp,nchp,1,nn)        call pastit (swet3,swetdeep,istrip,nchp,nchptot,1,nn)
1064        call pastit (capacity,capac,istrip,nchp,nchp,1,nn)        call pastit (capacity,capac,istrip,nchp,nchptot,1,nn)
1065        call pastit (snowdepth,snodep,istrip,nchp,nchp,1,nn)        call pastit (snowdepth,snodep,istrip,nchp,nchptot,1,nn)
1066    
1067  c**********************************************************************  c**********************************************************************
1068  c  Now update the theta and sh profiles with the new ground temperature  c  Now update the theta and sh profiles with the new ground temperature
# Line 966  c*************************************** Line 1077  c***************************************
1077        enddo        enddo
1078        enddo        enddo
1079    
1080          if(iudiag4.gt.0) then
1081          print *,' Bumping UDIAG4 '
1082          call paste2grd(th,igrd,chfrstr,istrip,nchp,
1083         1                   tempor4,ijall,nlay,nn,.false.)
1084          endif
1085    
1086    
1087        do i =1,istrip        do i =1,istrip
1088        sh(i,nlay+1) = qa(i)        sh(i,nlay+1) = qa(i)
1089        enddo        enddo
# Line 982  c*************************************** Line 1100  c***************************************
1100        enddo        enddo
1101        enddo        enddo
1102    
1103        if(tprof)then        if(2.eq.1)then
1104         CALL PNTPRF (1,IJALL,nlay,NYMD,NHMS,transth,'TRB T FLUX      ')        print *,' In turbio, just after updating th and sh - strip ',nn
1105         CALL PNTPRF (1,IJALL,nlay,NYMD,NHMS,transsh,'TRB Q FLUX      ')        print *,' field: th ',th
1106          print *,' field: sh ',sh
1107        endif        endif
1108    
1109  c tendency updates  c tendency updates
# Line 995  c ---------------- Line 1114  c ----------------
1114        do i =1,istrip        do i =1,istrip
1115        tends(i) = ( u(i,l)-tmpstrip(i) )        tends(i) = ( u(i,l)-tmpstrip(i) )
1116        enddo        enddo
1117        call pastit (tends,dumsc(1,l),istrip,nchp,nchp,1,nn)        call pastit (tends,dumsc(1,l),istrip,nchp,nchptot,1,nn)
1118    
1119        call strip2tile(vz(1,1,l),igrd,tmpstrip,nchp,ijall,        call strip2tile(vz(1,1,l),igrd,tmpstrip,nchp,ijall,
1120       1                                                 istrip,1,nn)       1                                                 istrip,1,nn)
1121        do i =1,istrip        do i =1,istrip
1122        tends(i) = ( v(i,l)-tmpstrip(i) )        tends(i) = ( v(i,l)-tmpstrip(i) )
1123        enddo        enddo
1124        call pastit (tends,dvmsc(1,l),istrip,nchp,nchp,1,nn)        call pastit (tends,dvmsc(1,l),istrip,nchp,nchptot,1,nn)
1125    
1126        call strip2tile(tz(1,1,l),igrd,tmpstrip,nchp,ijall,        call strip2tile(tz(1,1,l),igrd,tmpstrip,nchp,ijall,
1127       1                                                 istrip,1,nn)       1                                                 istrip,1,nn)
1128        do i =1,istrip        do i =1,istrip
1129        tends(i) = ( th(i,l)-tmpstrip(i) )        tends(i) = ( th(i,l)-tmpstrip(i) )
1130        enddo        enddo
1131        call pastit (tends,dtmsc(1,l),istrip,nchp,nchp,1,nn)  
1132          if(2.eq.1)then
1133          print *,' In turbio, tendencies for strip ',nn,' level ',l
1134          print *,' field: th ',tends
1135          endif
1136    
1137          call pastit (tends,dtmsc(1,l),istrip,nchp,nchptot,1,nn)
1138    
1139        call strip2tile(qz(1,1,l,1),igrd,tmpstrip,nchp,ijall,        call strip2tile(qz(1,1,l,1),igrd,tmpstrip,nchp,ijall,
1140       1                                                 istrip,1,nn)       1                                                 istrip,1,nn)
1141        do i =1,istrip        do i =1,istrip
1142        tends(i) = ( sh(i,l)-tmpstrip(i) )        tends(i) = ( sh(i,l)-tmpstrip(i) )
1143        enddo        enddo
       call pastit (tends,dqmsc(1,l,1),istrip,nchp,nchp,1,nn)  
1144    
1145        do nt = 1,ntracers-ptracers        if(2.eq.1)then
1146         call strip2tile(qz(1,1,L,ptracers+nt),igrd,tmpstrip,nchp,        print *,' In turbio, tendencies for strip ',nn,' level ',l
1147       1                                             ijall,istrip,1,nn)        print *,' field: sh ',tends
1148         do i =1,istrip        endif
1149          tends(i) = ( tracers(i,L,nt)-tmpstrip(i) )  
1150         enddo        call pastit (tends,dqmsc(1,l,1),istrip,nchp,nchptot,1,nn)
1151         call pastit (tends,dqmsc(1,L,ptracers+nt),istrip,nchp,nchp,1,nn)  
1152        enddo  c     do nt = 1,ntracers-ptracers
1153    c      call strip2tile(qz(1,1,L,ptracers+nt),igrd,tmpstrip,nchp,
1154    c    1                                             ijall,istrip,1,nn)
1155    c      do i =1,istrip
1156    c       tends(i) = ( tracers(i,L,nt)-tmpstrip(i) )
1157    c      enddo
1158    c      call pastit(tends,dqmsc(1,L,ptracers+nt),istrip,nchp,
1159    c    .                                     nchptot,1,nn)
1160    c     enddo
1161    
1162        enddo        enddo
1163    
# Line 1037  c note: the order, logic, and scaling of Line 1169  c note: the order, logic, and scaling of
1169  c       diagnostics is critical!  c       diagnostics is critical!
1170  c ------------------------------  c ------------------------------
1171    
1172        if(ievap.gt.0) then  c     if(ievap.gt.0) then
1173        do i=1,istrip  c     do i=1,istrip
1174        tmpstrip(i) = stqflux(i,nlay) * 86400  c     tmpstrip(i) = stqflux(i,nlay) * 86400
1175        enddo  c     enddo
1176        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1177       1                        qdiag(1,1,ievap),ijall,1,nn,.false.)  c    1                    qdiag(1,1,ievap,bi,bj),ijall,1,nn,.false.)
1178        endif  c     endif
1179        if(ieflux.gt.0) then  c     if(ieflux.gt.0) then
1180        do i=1,istrip  c     do i=1,istrip
1181        tmpstrip(i) = stqflux(i,nlay) * alhl  c     tmpstrip(i) = stqflux(i,nlay) * alhl
1182        enddo  c     enddo
1183        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1184       1                        qdiag(1,1,ieflux),ijall,1,nn,.false.)  c    1                    qdiag(1,1,ieflux,bi,bj),ijall,1,nn,.false.)
1185        endif  c     endif
1186        if(ihflux.gt.0) then  c     if(ihflux.gt.0) then
1187        do i=1,istrip  c     do i=1,istrip
1188        tmpstrip(i) = sttflux(i,nlay) * cp  c     tmpstrip(i) = sttflux(i,nlay) * cp
1189        enddo  c     enddo
1190        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1191       1                        qdiag(1,1,ihflux),ijall,1,nn,.false.)  c    1                    qdiag(1,1,ihflux,bi,bj),ijall,1,nn,.false.)
1192        endif  c     endif
1193        if(ituflux.gt.0) then  c     if(ituflux.gt.0) then
1194        call paste2grd(stuflux,igrd,chfrstr,istrip,nchp,  c     call paste2grd(stuflux,igrd,chfrstr,istrip,nchp,
1195       1                       qdiag(1,1,ituflux),ijall,nlay,nn,.false.)  c    1                   qdiag(1,1,ituflux,bi,bj),ijall,nlay,nn,.false.)
1196        endif  c     endif
1197        if(itvflux.gt.0) then  c     if(itvflux.gt.0) then
1198        call paste2grd(stvflux,igrd,chfrstr,istrip,nchp,  c     call paste2grd(stvflux,igrd,chfrstr,istrip,nchp,
1199       1                       qdiag(1,1,itvflux),ijall,nlay,nn,.false.)  c    1                   qdiag(1,1,itvflux,bi,bj),ijall,nlay,nn,.false.)
1200        endif  c     endif
1201        if(ittflux.gt.0) then  c     if(ittflux.gt.0) then
1202        do l=1,nlay  c     do l=1,nlay
1203        do i=1,istrip  c     do i=1,istrip
1204        sttflux(i,l) = sttflux(i,l) * cp  c     sttflux(i,l) = sttflux(i,l) * cp
1205        enddo  c     enddo
1206        enddo  c     enddo
1207        call paste2grd(sttflux,igrd,chfrstr,istrip,nchp,  c     call paste2grd(sttflux,igrd,chfrstr,istrip,nchp,
1208       1                       qdiag(1,1,ittflux),ijall,nlay,nn,.false.)  c    1                   qdiag(1,1,ittflux,bi,bj),ijall,nlay,nn,.false.)
1209        endif  c     endif
1210        if(itqflux.gt.0) then  c     if(itqflux.gt.0) then
1211        do l=1,nlay  c     do l=1,nlay
1212        do i=1,istrip  c     do i=1,istrip
1213        stqflux(i,l) = stqflux(i,l) * alhl  c     stqflux(i,l) = stqflux(i,l) * alhl
1214        enddo  c     enddo
1215        enddo  c     enddo
1216        call paste2grd(stqflux,igrd,chfrstr,istrip,nchp,  c     call paste2grd(stqflux,igrd,chfrstr,istrip,nchp,
1217       1                       qdiag(1,1,itqflux),ijall,nlay,nn,.false.)  c    1                   qdiag(1,1,itqflux,bi,bj),ijall,nlay,nn,.false.)
1218        endif  c     endif
1219        if(iri.gt.0) call paste2grd(sri,igrd,chfrstr,istrip,nchp,  c     if(iri.gt.0) call paste2grd(sri,igrd,chfrstr,istrip,nchp,
1220       1                           qdiag(1,1,iri),ijall,nlay,nn,.false.)  c    1                       qdiag(1,1,iri,bi,bj),ijall,nlay,nn,.false.)
1221        if(ikh.gt.0) call paste2grd(skh,igrd,chfrstr,istrip,nchp,  c     if(ikh.gt.0) call paste2grd(skh,igrd,chfrstr,istrip,nchp,
1222       1                           qdiag(1,1,ikh),ijall,nlay,nn,.false.)  c    1                       qdiag(1,1,ikh,bi,bj),ijall,nlay,nn,.false.)
1223        if(ikm.gt.0) call paste2grd(skm,igrd,chfrstr,istrip,nchp,  c     if(ikm.gt.0) call paste2grd(skm,igrd,chfrstr,istrip,nchp,
1224       1                           qdiag(1,1,ikm),ijall,nlay,nn,.false.)  c    1                       qdiag(1,1,ikm,bi,bj),ijall,nlay,nn,.false.)
1225        if(ict.gt.0) then  c     if(ict.gt.0) then
1226        call paste2grd(sct,igrd,chfrstr,istrip,nchp,  c     call paste2grd(sct,igrd,chfrstr,istrip,nchp,
1227       1                   qdiag(1,1,ict),ijall,1,nn,.false.)  c    1                   qdiag(1,1,ict,bi,bj),ijall,1,nn,.false.)
1228        endif  c     endif
1229        if(icu.gt.0) then  c     if(icu.gt.0) then
1230        call paste2grd(scu,igrd,chfrstr,istrip,nchp,  c     call paste2grd(scu,igrd,chfrstr,istrip,nchp,
1231       1                   qdiag(1,1,icu),ijall,1,nn,.false.)  c    1                   qdiag(1,1,icu,bi,bj),ijall,1,nn,.false.)
1232        endif  c     endif
1233        if(iwinds.gt.0) then  c     if(iwinds.gt.0) then
1234        call paste2grd(swinds,igrd,chfrstr,istrip,nchp,  c     call paste2grd(swinds,igrd,chfrstr,istrip,nchp,
1235       1                      qdiag(1,1,iwinds),ijall,1,nn,.false.)  c    1                      qdiag(1,1,iwinds,bi,bj),ijall,1,nn,.false.)
1236        endif  c     endif
1237        if(iuflux.gt.0) then  c     if(iuflux.gt.0) then
1238        call paste2grd(stuflux(1,nlay),igrd,chfrstr,istrip,nchp,  c     call paste2grd(stuflux(1,nlay),igrd,chfrstr,istrip,nchp,
1239       1                       qdiag(1,1,iuflux),ijall,1,nn,.false.)  c    1                       qdiag(1,1,iuflux,bi,bj),ijall,1,nn,.false.)
1240        endif  c     endif
1241        if(ivflux.gt.0) then  c     if(ivflux.gt.0) then
1242        call paste2grd(stvflux(1,nlay),igrd,chfrstr,istrip,nchp,  c     call paste2grd(stvflux(1,nlay),igrd,chfrstr,istrip,nchp,
1243       1                       qdiag(1,1,ivflux),ijall,1,nn,.false.)  c    1                       qdiag(1,1,ivflux,bi,bj),ijall,1,nn,.false.)
1244        endif  c     endif
1245        if(iustar.gt.0) then  c     if(iustar.gt.0) then
1246        call paste2grd(sustar,igrd,chfrstr,istrip,nchp,  c     call paste2grd(sustar,igrd,chfrstr,istrip,nchp,
1247       1                      qdiag(1,1,iustar),ijall,1,nn,.false.)  c    1                      qdiag(1,1,iustar,bi,bj),ijall,1,nn,.false.)
1248        endif  c     endif
1249        if(iz0.gt.0) then  c     if(iz0.gt.0) then
1250        call paste2grd(sz0,igrd,chfrstr,istrip,nchp,  c     call paste2grd(sz0,igrd,chfrstr,istrip,nchp,
1251       1                   qdiag(1,1,iz0),ijall,1,nn,.false.)  c    1                   qdiag(1,1,iz0,bi,bj),ijall,1,nn,.false.)
1252        endif  c     endif
1253        if(ifrqtrb.gt.0) then  c     if(ifrqtrb.gt.0) then
1254        call paste2grd(frqtrb,igrd,chfrstr,istrip,nchp,  c     call paste2grd(frqtrb,igrd,chfrstr,istrip,nchp,
1255       1                      qdiag(1,1,ifrqtrb),ijall,1,nn,.false.)  c    1                      qdiag(1,1,ifrqtrb,bi,bj),ijall,1,nn,.false.)
1256        endif  c     endif
1257        if(ipbl.gt.0) then  c     if(ipbl.gt.0) then
1258        call paste2grd(pbldpth,igrd,chfrstr,istrip,nchp,  c     call paste2grd(pbldpth,igrd,chfrstr,istrip,nchp,
1259       1                       qdiag(1,1,ipbl),ijall,1,nn,.false.)  c    1                       qdiag(1,1,ipbl,bi,bj),ijall,1,nn,.false.)
1260        endif  c     endif
1261        if(iu2m.gt.0) then  c     if(iu2m.gt.0) then
1262        call paste2grd(stu2m,igrd,chfrstr,istrip,nchp,  c     call paste2grd(stu2m,igrd,chfrstr,istrip,nchp,
1263       1                     qdiag(1,1,iu2m),ijall,1,nn,.true.)  c    1                     qdiag(1,1,iu2m,bi,bj),ijall,1,nn,.true.)
1264        endif  c     endif
1265        if(iv2m.gt.0) then  c     if(iv2m.gt.0) then
1266        call paste2grd(stv2m,igrd,chfrstr,istrip,nchp,  c     call paste2grd(stv2m,igrd,chfrstr,istrip,nchp,
1267       1                     qdiag(1,1,iv2m),ijall,1,nn,.true.)  c    1                     qdiag(1,1,iv2m,bi,bj),ijall,1,nn,.true.)
1268        endif  c     endif
1269        if(it2m.gt.0) then  c     if(it2m.gt.0) then
1270        call paste2grd(stt2m,igrd,chfrstr,istrip,nchp,  c     call paste2grd(stt2m,igrd,chfrstr,istrip,nchp,
1271       1                     qdiag(1,1,it2m),ijall,1,nn,.true.)  c    1                     qdiag(1,1,it2m,bi,bj),ijall,1,nn,.true.)
1272        endif  c     endif
1273        if(iq2m.gt.0) then  c     if(iq2m.gt.0) then
1274        do i=1,istrip  c     do i=1,istrip
1275           if( stq2m(i).ne.undef ) then  c        if( stq2m(i).ne.undef ) then
1276               tmpstrip(i) = stq2m(i) * 1000  c            tmpstrip(i) = stq2m(i) * 1000
1277           else  c        else
1278               tmpstrip(i) = undef  c            tmpstrip(i) = undef
1279           endif  c        endif
1280        enddo  c     enddo
1281        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1282       1                     qdiag(1,1,iq2m),ijall,1,nn,.true.)  c    1                     qdiag(1,1,iq2m,bi,bj),ijall,1,nn,.true.)
1283        endif  c     endif
1284        if(iu10m.gt.0) then  c     if(iu10m.gt.0) then
1285        call paste2grd(stu10m,igrd,chfrstr,istrip,nchp,  c     call paste2grd(stu10m,igrd,chfrstr,istrip,nchp,
1286       1                      qdiag(1,1,iu10m),ijall,1,nn,.false.)  c    1                      qdiag(1,1,iu10m,bi,bj),ijall,1,nn,.false.)
1287        endif  c     endif
1288        if(iv10m.gt.0) then  c     if(iv10m.gt.0) then
1289        call paste2grd(stv10m,igrd,chfrstr,istrip,nchp,  c     call paste2grd(stv10m,igrd,chfrstr,istrip,nchp,
1290       1                      qdiag(1,1,iv10m),ijall,1,nn,.false.)  c    1                      qdiag(1,1,iv10m,bi,bj),ijall,1,nn,.false.)
1291        endif  c     endif
1292        if(it10m.gt.0) then  c     if(it10m.gt.0) then
1293        call paste2grd(stt10m,igrd,chfrstr,istrip,nchp,  c     call paste2grd(stt10m,igrd,chfrstr,istrip,nchp,
1294       1                      qdiag(1,1,it10m),ijall,1,nn,.false.)  c    1                      qdiag(1,1,it10m,bi,bj),ijall,1,nn,.false.)
1295        endif  c     endif
1296        if(iq10m.gt.0) then  c     if(iq10m.gt.0) then
1297        do i=1,istrip  c     do i=1,istrip
1298           if( stq10m(i).ne.undef ) then  c        if( stq10m(i).ne.undef ) then
1299               tmpstrip(i) = stq10m(i) * 1000  c            tmpstrip(i) = stq10m(i) * 1000
1300           else  c        else
1301               tmpstrip(i) = undef  c            tmpstrip(i) = undef
1302           endif  c        endif
1303        enddo  c     enddo
1304        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1305       1                      qdiag(1,1,iq10m),ijall,1,nn,.false.)  c    1                      qdiag(1,1,iq10m,bi,bj),ijall,1,nn,.false.)
1306        endif  c     endif
   
 c Fill 2m and 10m Couplings (Not Re-producible)  
 c ---------------------------------------------  
       call paste2grd (  stu2m,igrd,chfrstr,istrip,nchp, u2m,ijall,1,  
      .                                                      nn,.true. )  
       call paste2grd (  stv2m,igrd,chfrstr,istrip,nchp, v2m,ijall,1,  
      .                                                      nn,.true. )  
       call paste2grd (  stt2m,igrd,chfrstr,istrip,nchp, t2m,ijall,1,  
      .                                                      nn,.true. )  
       call paste2grd (  stq2m,igrd,chfrstr,istrip,nchp, q2m,ijall,1,  
      .                                                      nn,.true. )  
       call paste2grd ( stu10m,igrd,chfrstr,istrip,nchp,u10m,ijall,1,  
      .                                                      nn,.false.)  
       call paste2grd ( stv10m,igrd,chfrstr,istrip,nchp,v10m,ijall,1,  
      .                                                      nn,.false.)  
       call paste2grd ( stt10m,igrd,chfrstr,istrip,nchp,t10m,ijall,1,  
      .                                                      nn,.false.)  
       call paste2grd ( stq10m,igrd,chfrstr,istrip,nchp,q10m,ijall,1,  
      .                                                      nn,.false.)  
1307    
1308  c**********************************************************************  c**********************************************************************
1309  c   more diagnostics: land surface model parameters  c   more diagnostics: land surface model parameters
1310  c**********************************************************************  c**********************************************************************
1311    
1312        if(itdeep.gt.0)call paste2grd(td,igrd,chfrstr,istrip,nchp,  c     if(itdeep.gt.0)call paste2grd(td,igrd,chfrstr,istrip,nchp,
1313       .                         qdiag(1,1,itdeep),ijall,1,nn,.false.)  c    .                     qdiag(1,1,itdeep,bi,bj),ijall,1,nn,.false.)
1314        if(iqcanopy .gt.0)call paste2grd(qa,igrd,chfrstr,istrip,nchp,  c     if(iqcanopy .gt.0)call paste2grd(qa,igrd,chfrstr,istrip,nchp,
1315       .                      qdiag(1,1,iqcanopy) ,ijall,1,nn,.false.)  c    .                  qdiag(1,1,iqcanopy,bi,bj) ,ijall,1,nn,.false.)
1316        if(ismshal  .gt.0)call paste2grd(swet1,igrd,chfrstr,istrip,nchp,  c     if(ismshal  .gt.0)call paste2grd(swet1,igrd,chfrstr,istrip,nchp,
1317       .                      qdiag(1,1,ismshal)  ,ijall,1,nn,.false.)  c    .                  qdiag(1,1,ismshal,bi,bj)  ,ijall,1,nn,.false.)
1318        if(ismroot  .gt.0)call paste2grd(swet2,igrd,chfrstr,istrip,nchp,  c     if(ismroot  .gt.0)call paste2grd(swet2,igrd,chfrstr,istrip,nchp,
1319       .                      qdiag(1,1,ismroot)  ,ijall,1,nn,.false.)  c    .                  qdiag(1,1,ismroot,bi,bj)  ,ijall,1,nn,.false.)
1320        if(ismdeep  .gt.0)call paste2grd(swet3,igrd,chfrstr,istrip,nchp,  c     if(ismdeep  .gt.0)call paste2grd(swet3,igrd,chfrstr,istrip,nchp,
1321       .                      qdiag(1,1,ismdeep)  ,ijall,1,nn,.false.)  c    .                  qdiag(1,1,ismdeep,bi,bj)  ,ijall,1,nn,.false.)
1322        if(icapacity.gt.0)call paste2grd(capacity,igrd,chfrstr,istrip,  c     if(icapacity.gt.0)call paste2grd(capacity,igrd,chfrstr,istrip,
1323       .                 nchp,qdiag(1,1,icapacity),ijall,1,nn,.false.)  c    .             nchp,qdiag(1,1,icapacity,bi,bj),ijall,1,nn,.false.)
1324        if(isnow.gt.0)call paste2grd(snowdepth,igrd,chfrstr,istrip,nchp,  c     if(isnow.gt.0)call paste2grd(snowdepth,igrd,chfrstr,istrip,nchp,
1325       .                      qdiag(1,1,isnow)    ,ijall,1,nn,.false.)  c    .                  qdiag(1,1,isnow,bi,bj)    ,ijall,1,nn,.false.)
   
 c**********************************************************************  
       IF(Iudiag1.GT.0) then  
       call paste2grd(checktrb,igrd,chfrstr,istrip,nchp,  
      1                       qdiag(1,1,iudiag1),ijall,nlay,nn,.false.)  
       endif  
1326    
1327  c**********************************************************************  c**********************************************************************
1328  c end regions loop  c end regions loop
# Line 1231  c -------------------------------------- Line 1338  c --------------------------------------
1338    
1339  c prevent ice or snow from melting  c prevent ice or snow from melting
1340  c ---------------------------------------------------------------------  c ---------------------------------------------------------------------
1341        do i =1,nchp        do i =1,nchptot
1342        if( (icethk(i).gt.0.).and.(tground(i).gt.tice) ) tground(i)=tice        if( (icethk(i).gt.0.).and.(tground(i).gt.tice) ) tground(i)=tice
1343        enddo        enddo
1344    
# Line 1264  c ----------------------------------- Line 1371  c -----------------------------------
1371  C Return Tendencies and Couplings to Grid Space  C Return Tendencies and Couplings to Grid Space
1372  c ---------------------------------------------  c ---------------------------------------------
1373        do l = 1,nlay        do l = 1,nlay
1374        call msc2grd(igrd,chfr,dumsc(1,L),nchp,nchp,fracland,        call msc2grd(igrd,chfr,dumsc(1,L),nchp,nchptot,fracland,
1375       .                                              duturb(1,1,L),im,jm)       .                                              duturb(1,1,L),im,jm)
1376        call msc2grd(igrd,chfr,dvmsc(1,L),nchp,nchp,fracland,        call msc2grd(igrd,chfr,dvmsc(1,L),nchp,nchptot,fracland,
1377       .                                              dvturb(1,1,L),im,jm)       .                                              dvturb(1,1,L),im,jm)
1378        call msc2grd(igrd,chfr,dtmsc(1,L),nchp,nchp,fracland,        call msc2grd(igrd,chfr,dtmsc(1,L),nchp,nchptot,fracland,
1379       .                                              dtturb(1,1,L),im,jm)       .                                              dtturb(1,1,L),im,jm)
1380        do nt = 1,ntracers        do nt = 1,ntracers
1381        call msc2grd(igrd,chfr,dqmsc(1,L,nt),nchp,nchp,fracland,        call msc2grd(igrd,chfr,dqmsc(1,L,nt),nchp,nchptot,fracland,
1382       .                                           dqturb(1,1,L,nt),im,jm)       .                                           dqturb(1,1,L,nt),im,jm)
1383        enddo        enddo
1384        call msc2grd(igrd,chfr,  tke(1,L),nchp,nchp,fracland,        call msc2grd(igrd,chfr,  tke(1,L),nchp,nchptot,fracland,
1385       .                                              qqgrid(1,1,L),im,jm)       .                                              qqgrid(1,1,L),im,jm)
1386    
1387        call msc2grd(igrd,chfr, fccmsc(1,L),nchp,nchp,fracland,        call msc2grd(igrd,chfr, fccmsc(1,L),nchp,nchptot,fracland,
1388       .                                              fcctmp(1,1,L),im,jm)       .                                              fcctmp(1,1,L),im,jm)
1389        call msc2grd(igrd,chfr,qliqmsc(1,L),nchp,nchp,fracland,        call msc2grd(igrd,chfr,qliqmsc(1,L),nchp,nchptot,fracland,
1390       .                                             qliqtmp(1,1,L),im,jm)       .                                             qliqtmp(1,1,L),im,jm)
1391        enddo        enddo
1392    
# Line 1302  c -------------------------------------- Line 1409  c --------------------------------------
1409         if (itrbfcc.gt.0) then         if (itrbfcc.gt.0) then
1410         do j=1,jm         do j=1,jm
1411         do i=1,im         do i=1,im
1412         qdiag(i,j,itrbfcc+L-1) =  qdiag(i,j,itrbfcc+L-1) + fcctmp(i,j,L)         qdiag(i,j,itrbfcc+L-1,bi,bj) =  qdiag(i,j,itrbfcc+L-1,bi,bj) +
1413         .                                              fcctmp(i,j,L)
1414         enddo         enddo
1415         enddo         enddo
1416         endif         endif
# Line 1310  c -------------------------------------- Line 1418  c --------------------------------------
1418         if (itrbqliq.gt.0) then         if (itrbqliq.gt.0) then
1419         do j=1,jm         do j=1,jm
1420         do i=1,im         do i=1,im
1421         qdiag(i,j,itrbqliq+L-1)=qdiag(i,j,itrbqliq+L-1)+         qdiag(i,j,itrbqliq+L-1,bi,bj)=qdiag(i,j,itrbqliq+L-1,bi,bj)+
1422       .                                             qliqtmp(i,j,L)*1.e6       .                                             qliqtmp(i,j,L)*1.e6
1423         enddo         enddo
1424         enddo         enddo
# Line 1323  C Ground Temperature, snow depth and sha Line 1431  C Ground Temperature, snow depth and sha
1431        do j = 1,jm        do j = 1,jm
1432         do i = 1,im         do i = 1,im
1433           tgz(i,j) = 0.           tgz(i,j) = 0.
         snow(i,j) = 0.  
         gwet(i,j) = 0.  
1434         enddo         enddo
1435        enddo        enddo
1436        call msc2grd(igrd,chfr,tground ,nchp,nchp,fracland,tgz ,im,jm)        call msc2grd(igrd,chfr,tground ,nchp,nchptot,fracland,tgz ,im,jm)
       call msc2grd(igrd,chfr,snodep  ,nchp,nchp,fracland,snow,im,jm)  
       call msc2grd(igrd,chfr,swetshal,nchp,nchp,fracland,gwet,im,jm)  
1437    
1438  c *********************************************************************  c *********************************************************************
1439  c **** increment diagnostic array for ground and surface temperatures,  c **** increment diagnostic array for ground and surface temperatures,
1440  c ***       ground temp tendency, and ground humidity  c ***       ground temp tendency, and ground humidity
1441  c *********************************************************************  c *********************************************************************
1442    
1443        if(itground.gt.0) then  c     if(itground.gt.0) then
1444        do i =1,ijall  c     do j =1,jm
1445        qdiag(i,1,itground) = qdiag(i,1,itground) + tgz(i,1)  c     do i =1,im
1446        enddo  c     qdiag(i,j,itground,bi,bj) = qdiag(i,j,itground,bi,bj) + tgz(i,j)
1447        endif  c     enddo
1448    c     enddo
1449    c     endif
1450    
1451        if(itcanopy.gt.0) then        if(itcanopy.gt.0) then
1452        do i =1,ijall        do j =1,jm
1453        qdiag(i,1,itcanopy) = qdiag(i,1,itcanopy) + tgz(i,1)        do i =1,im
1454          qdiag(i,j,itcanopy,bi,bj) = qdiag(i,j,itcanopy,bi,bj) + tgz(i,j)
1455          enddo
1456        enddo        enddo
1457        endif        endif
1458    
1459        if(its.gt.0) then        if(its.gt.0) then
1460        do i =1,ijall        do j =1,jm
1461        tmpstrip(i) = tz(i,1,nlay) * pkht(i,1,nlay)        do i =1,im
1462          tmpdiag(i,j) = tz(i,j,nlay) * pkht(i,j,nlay)
1463          enddo
1464          enddo
1465          do j =1,jm
1466          do i =1,im
1467          qdiag(i,j,its,bi,bj) = qdiag(i,j,its,bi,bj) + tmpdiag(i,j)
1468        enddo        enddo
       do i =1,ijall  
       qdiag(i,1,its) = qdiag(i,1,its) + tmpstrip(i)  
1469        enddo        enddo
1470        endif        endif
1471    
1472        if(idtg.gt.0) then        if(idtg.gt.0) then
1473        do i =1,ijall        do j =1,jm
1474        qdiag(i,1,idtg) = qdiag(i,1,idtg) + tgz(i,1)        do i =1,im
1475          qdiag(i,j,idtg,bi,bj) = qdiag(i,j,idtg,bi,bj) + tgz(i,j)
1476          enddo
1477        enddo        enddo
1478        endif        endif
1479    
# Line 1369  c ************************************** Line 1483  c **************************************
1483        do L = 1,nlay        do L = 1,nlay
1484    
1485         if(iturbu.gt.0) then         if(iturbu.gt.0) then
1486         do i =1,ijall         do j =1,jm
1487          qdiag(i,1,iturbu+l-1) =  qdiag(i,1,iturbu+l-1)         do i =1,im
1488       .                      + duturb(i,1,l) * atimstp * secday          qdiag(i,j,iturbu+l-1,bi,bj) =  qdiag(i,j,iturbu+l-1,bi,bj)
1489         .                      + duturb(i,j,l) * atimstp * secday
1490           enddo
1491         enddo         enddo
1492         endif         endif
1493    
1494         if(iturbv.gt.0) then         if(iturbv.gt.0) then
1495         do i =1,ijall         do j =1,jm
1496          qdiag(i,1,iturbv+l-1) =  qdiag(i,1,iturbv+l-1)         do i =1,im
1497       .                      + dvturb(i,1,l) * atimstp * secday          qdiag(i,j,iturbv+l-1,bi,bj) =  qdiag(i,j,iturbv+l-1,bi,bj)
1498         .                      + dvturb(i,j,l) * atimstp * secday
1499           enddo
1500         enddo         enddo
1501         endif         endif
1502    
1503         if(iturbq.gt.0) then         if(iturbq.gt.0) then
1504         do i =1,ijall         do j =1,jm
1505          qdiag(i,1,iturbq+l-1) =  qdiag(i,1,iturbq+l-1)         do i =1,im
1506       .                      + dqturb(i,1,l,1) * atimstp * secday * 1000          qdiag(i,j,iturbq+l-1,bi,bj) =  qdiag(i,j,iturbq+l-1,bi,bj)
1507         .                      + dqturb(i,j,l,1) * atimstp * secday * 1000
1508           enddo
1509         enddo         enddo
1510         endif         endif
1511    
1512         if(iturbt.gt.0) then         if(iturbt.gt.0) then
1513         do i =1,ijall         do j =1,jm
1514          qdiag(i,1,iturbt+l-1) =  qdiag(i,1,iturbt+l-1)         do i =1,im
1515       .                   + dtturb(i,1,l) * pkz(i,1,l)*atimstp*secday          qdiag(i,j,iturbt+l-1,bi,bj) =  qdiag(i,j,iturbt+l-1,bi,bj)
1516         .                   + dtturb(i,j,l) * pkz(i,1,l)*atimstp*secday
1517           enddo
1518           enddo
1519           endif
1520    
1521           if(iudiag1.gt.0) then
1522           do j =1,jm
1523           do i =1,im
1524            qdiag(i,j,iudiag1+l-1,bi,bj) =  qdiag(i,j,iudiag1+l-1,bi,bj)
1525         .                      + tempor1(i,j,l)
1526           enddo
1527           enddo
1528           endif
1529    
1530           if(iudiag2.gt.0) then
1531           do j =1,jm
1532           do i =1,im
1533            qdiag(i,j,iudiag2+l-1,bi,bj) =  qdiag(i,j,iudiag2+l-1,bi,bj)
1534         .                      + tempor2(i,j,l)
1535           enddo
1536           enddo
1537           endif
1538    
1539           if(iudiag3.gt.0) then
1540           do j =1,jm
1541           do i =1,im
1542            qdiag(i,j,iudiag3+l-1,bi,bj) =  qdiag(i,j,iudiag3+l-1,bi,bj)
1543         .                      + tempor3(i,j,l)
1544           enddo
1545           enddo
1546           endif
1547    
1548           if(iudiag4.gt.0) then
1549           do j =1,jm
1550           do i =1,im
1551            qdiag(i,j,iudiag4+l-1,bi,bj) =  qdiag(i,j,iudiag4+l-1,bi,bj)
1552         .                      + tempor4(i,j,l)
1553           enddo
1554         enddo         enddo
1555         endif         endif
1556    
1557        enddo        enddo
1558    
1559           if(isdiag1.gt.0) then
1560           do j =1,jm
1561           do i =1,im
1562            qdiag(i,j,isdiag1,bi,bj) =  qdiag(i,j,isdiag1,bi,bj)
1563         .                      + tempsfc1(i,j)
1564           enddo
1565           enddo
1566           endif
1567    
1568           if(isdiag2.gt.0) then
1569           do j =1,jm
1570           do i =1,im
1571            qdiag(i,j,isdiag2,bi,bj) =  qdiag(i,j,isdiag2,bi,bj)
1572         .                      + tempsfc2(i,j)
1573           enddo
1574           enddo
1575           endif
1576    
1577  c pi-weight the theta and moisture tendencies  c pi-weight the theta and moisture tendencies
1578  c -------------------------------------------  c -------------------------------------------
1579        do i =1,ijall        do i =1,ijall
# Line 1407  c -------------------------------------- Line 1583  c --------------------------------------
1583         do i =1,ijall         do i =1,ijall
1584          duturb(i,1,l) = duturb(i,1,l)*atimstp          duturb(i,1,l) = duturb(i,1,l)*atimstp
1585          dvturb(i,1,l) = dvturb(i,1,l)*atimstp          dvturb(i,1,l) = dvturb(i,1,l)*atimstp
1586          dtturb(i,1,l) = dtturb(i,1,l)*thtgz(i)  c       dtturb(i,1,l) = dtturb(i,1,l)*thtgz(i)
1587            dtturb(i,1,l) = dtturb(i,1,l)*atimstp
1588         enddo         enddo
1589         do nt = 1,ntracers         do nt = 1,ntracers
1590          do i =1,ijall          do i =1,ijall
1591           dqturb(i,1,l,nt) = dqturb(i,1,l,nt)*thtgz(i)  c        dqturb(i,1,l,nt) = dqturb(i,1,l,nt)*thtgz(i)
1592             dqturb(i,1,l,nt) = dqturb(i,1,l,nt)*atimstp
1593          enddo          enddo
1594         enddo         enddo
1595        enddo        enddo
# Line 1434  c ************************************** Line 1612  c **************************************
1612  c ****                bump diagnostic counters                      ***  c ****                bump diagnostic counters                      ***
1613  c *********************************************************************  c *********************************************************************
1614    
1615    #ifdef ALLOW_DIAGNOSTICS
1616          if( (bi.eq.1) .and. (bj.eq.1) ) then
1617        nturbt    = nturbt    + 1        nturbt    = nturbt    + 1
1618        nturbq    = nturbq    + 1        nturbq    = nturbq    + 1
1619        nturbu    = nturbu    + 1        nturbu    = nturbu    + 1
# Line 1512  c ************************************** Line 1692  c **************************************
1692        ntrbqliq  = ntrbqliq  + 1        ntrbqliq  = ntrbqliq  + 1
1693        ntrbfcc   = ntrbfcc   + 1        ntrbfcc   = ntrbfcc   + 1
1694    
1695          nsdiag1 = nsdiag1 + 1
1696          nsdiag2 = nsdiag2 + 1
1697          nudiag1 = nudiag1 + 1
1698          nudiag2 = nudiag2 + 1
1699          nudiag3 = nudiag3 + 1
1700          nudiag4 = nudiag4 + 1
1701          endif
1702    
1703    #endif
1704    
1705        return        return
1706        end        end
1707        SUBROUTINE TRBFLX (NN,TH,THV,SH,U,V,QQ,PL,PLE,PLK,PLKE,DPSTR,        SUBROUTINE TRBFLX (NN,TH,THV,SH,U,V,QQ,PL,PLE,PLK,PLKE,DPSTR,
# Line 1547  C    Z0            -         SURFACE ROU Line 1737  C    Z0            -         SURFACE ROU
1737  C    tracers       -         array of passive tracers  C    tracers       -         array of passive tracers
1738  C    ntrace        -         number of tracers to be diffused  C    ntrace        -         number of tracers to be diffused
1739  C    ntracedim     -         outer dimension of tracers array  C    ntracedim     -         outer dimension of tracers array
1740  C    DTAU          -         TIME CHANGE PER ITERATION OF TRBLFX  C    DTAU          -         TIME CHANGE PER ITERATION OF TRBFLX
1741  C    ITRTRB        -         NUMBER OF ITERATIONS OF TRBLFX  C    ITRTRB        -         NUMBER OF ITERATIONS OF TRBFLX
1742  C    KMBG          -         BACKGROUND VALUE OF MOMENTUM TRANSFER COEF  C    KMBG          -         BACKGROUND VALUE OF MOMENTUM TRANSFER COEF
1743  C    KHBG          -         BACKGROUND VALUE OF HEAT TRANSFER COEF  C    KHBG          -         BACKGROUND VALUE OF HEAT TRANSFER COEF
1744  C    NLEV          -         NUMBER OF ATMOSPHERIC LEVELS TO CALCULATE  C    NLEV          -         NUMBER OF ATMOSPHERIC LEVELS TO CALCULATE
# Line 1560  C     ------- Line 1750  C     -------
1750  C    PROFILES RETURNED WITH UPDATED VALUES  C    PROFILES RETURNED WITH UPDATED VALUES
1751  C  C
1752  C**********************************************************************  C**********************************************************************
1753  C        implicit none
1754  C  
1755    C Argument list declarations
1756          integer nn,irun,nlev,ntrace,ntracedim,itrtrb,nhms,nymd
1757          _RL TH(irun,NLEV+1),THV(irun,NLEV+1),SH(irun,NLEV+1)
1758          _RL U(irun,NLEV+1),V(irun,NLEV+1),QQ(irun,NLEV)
1759          _RL PL(irun,NLEV),PLE(irun,NLEV+1),PLK(irun,NLEV)
1760          _RL PLKE(irun,NLEV+1),DPSTR(irun,NLEV)
1761          integer IWATER(irun)
1762          _RL Z0(irun)
1763          _RL tracers(irun,nlev+1,ntracedim)
1764          _RL dtau,KMBG,KHBG
1765          LOGICAL QBEG,TPROF
1766          _RL SWINDS(irun)
1767          _RL SRI(irun,nlev), ET(irun,nlev)
1768          _RL EU (irun,nlev)
1769          _RL WU(irun,nlev)
1770          _RL WV (irun,nlev), pbldpth(irun)
1771          _RL sustar(irun), sz0(irun)
1772          _RL freqdg(irun,nlev-1)
1773          _RL sct(irun), scu(irun)
1774          _RL stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)
1775          _RL stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)
1776          _RL grav,cp,rgas,faceps,virtcon,undef
1777          _RL eturb(irun),dedqa(irun),dedtc(irun)
1778          _RL hsturb(irun),dhsdqa(irun),dhsdtc(irun)
1779          _RL dshdthg(irun,nlev),dthdthg(irun,nlev)
1780          _RL dshdshg(irun,nlev),dthdshg(irun,nlev)
1781          _RL transth(irun,nlev),transsh(irun,nlev)
1782          _RL ctsave(irun),xxsave(irun),yysave(irun)
1783          _RL zetasave(irun),xlsave(irun,nlev),khsave(irun,nlev)
1784          _RL qliq(irun,nlev),turbfcc(irun,nlev)
1785    
1786    C Local Variables
1787          _RL b1,b3,alpha,halpha,qqmin,qbustr
1788        PARAMETER ( B1      =   16.6    )          PARAMETER ( B1      =   16.6    )  
1789        PARAMETER ( B3      = 1. / B1  )          PARAMETER ( B3      = 1. / B1  )  
1790        PARAMETER ( ALPHA   = 0.1       )        PARAMETER ( ALPHA   = 0.1       )
1791        PARAMETER ( HALPHA = ALPHA * 0.5 )        PARAMETER ( HALPHA = ALPHA * 0.5 )
1792        PARAMETER ( QQMIN  = 0.005      )        PARAMETER ( QQMIN  = 0.005      )
1793        PARAMETER ( QBUSTR = 2.550952   )        PARAMETER ( QBUSTR = 2.550952   )
1794        real    argmax, onethrd, z1pem25, b2        _RL    argmax, onethrd, z1pem25, b2, two
1795        PARAMETER (ARGMAX = 30.)        PARAMETER (ARGMAX = 30.)
1796        PARAMETER (ONETHRD = 1./3. )        PARAMETER (ONETHRD = 1./3. )
1797        PARAMETER (Z1PEM25 = 1.E-25)        PARAMETER (Z1PEM25 = 1.E-25)
1798        PARAMETER ( B2    =  10.1 )        PARAMETER ( B2    =  10.1 )
1799        PARAMETER ( two   =   2.0 )        PARAMETER ( two   =   2.0 )
1800    
1801        DIMENSION TH(irun,NLEV+1),THV(irun,NLEV+1),SH(irun,NLEV+1)        _RL AHS (irun), HS(irun)
1802        DIMENSION U(irun,NLEV+1),V(irun,NLEV+1),QQ(irun,NLEV)        _RL XX  (irun), YY(irun), CU(irun)
1803        DIMENSION PL(irun,NLEV),PLE(irun,NLEV+1),PLK(irun,NLEV)        _RL CT(irun),  USTAR(irun)
1804        DIMENSION PLKE(irun,NLEV+1),DPSTR(irun,NLEV)        _RL RIB(irun),   ZETA(irun), WS(irun)
1805        DIMENSION IWATER(irun),Z0(irun)        _RL DTHS(irun), DELTHS(irun)
1806        real tracers(irun,nlev+1,ntracedim)        _RL DTHL(irun), DELTHL(irun)
1807        real KMBG,KHBG        _RL RIBIN(irun),CUIN(irun)
1808        LOGICAL QBEG,TPROF        _RL CTIN(irun),ZETAIN(irun)
1809        real eturb(irun),dedqa(irun),dedtc(irun)        _RL USTARIN(irun),RHOSIN(irun),Z0IN(irun)
1810        real hsturb(irun),dhsdqa(irun),dhsdtc(irun)        _RL qqcolmin(irun),qqcolmax(irun),levpbl(irun)
1811        real dshdthg(irun,nlev),dthdthg(irun,nlev)  
1812        real dshdshg(irun,nlev),dthdshg(irun,nlev)        _RL ADZ1(irun,nlev), DZ1TMP(irun,nlev)
1813        real transth(irun,nlev),transsh(irun,nlev)        _RL DZ3(irun,nlev), TEMP(irun,nlev)
1814        real ctsave(irun),xxsave(irun),yysave(irun)        _RL DV(irun,nlev), DTHV(irun,nlev)
1815        real zetasave(irun),xlsave(irun,nlev),khsave(irun,nlev)        _RL DPK(irun,nlev), STRT(irun,nlev)
1816        real qliq(irun,nlev),turbfcc(irun,nlev)        _RL DW2(irun,nlev), RI(irun,nlev)
1817          _RL RHOZPK(irun,nlev), Q(irun,nlev)
1818  C Diagnostic Variables        _RL RIINIT(irun,nlev), DU(irun,nlev)
1819  C --------------------        _RL QQINIT(irun,nlev), RHOKDZ(irun,nlev)
1820        DIMENSION SWINDS(irun)        _RL RHODZ2(irun,nlev)
1821        DIMENSION    SRI(irun,nlev), ET(irun,nlev)        _RL KM(irun,nlev), KH(irun,nlev)
1822        DIMENSION    EU (irun,nlev)  
1823        DIMENSION    WU(irun,nlev)        _RL DELTH  (irun,nlev+1), DELSH (irun,nlev+1)
1824        DIMENSION    WV (irun,nlev), pbldpth(irun)        _RL FLXFAC (irun,nlev+1)
1825        DIMENSION sustar(irun),          sz0(irun)        _RL FLXFPK (irun,nlev+1)
1826        DIMENSION    sct(irun),          scu(irun)  
1827        dimension stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)        _RL ADZ2   (irun,nlev-1), RHODZ1(irun,nlev-1)
1828        dimension stu10m(irun),stv10m(irun),        _RL VKZE   (irun,nlev-1), VKZM  (irun,nlev-1)
1829       1                              stt10m(irun),stq10m(irun)        _RL XL     (irun,nlev-1), QXLM  (irun,nlev-1)
1830        DIMENSION freqdg(irun,nlev-1)        _RL QQE    (irun,nlev-1), QE    (irun,nlev-1)
1831          _RL P3     (irun,nlev-1), XQ    (irun,nlev-1)
1832  C Dynamic Variables        _RL XLDIAG (irun,nlev-1), FLXFCE(irun,nlev-1)
 C -----------------  
       DIMENSION AHS (irun), HS(irun)  
       DIMENSION XX  (irun), YY(irun), CU(irun)  
       DIMENSION CT     (irun),  USTAR(irun)  
       DIMENSION RIB    (irun),   ZETA(irun), WS(irun)  
       DIMENSION DTHS   (irun), DELTHS(irun)  
       DIMENSION DTHL   (irun), DELTHL(irun)  
       DIMENSION TG     (irun)  
       DIMENSION RIBIN  (irun),   CUIN(irun)  
       DIMENSION CTIN   (irun), ZETAIN(irun)  
       DIMENSION USTARIN(irun), RHOSIN(irun), Z0IN(irun)  
       DIMENSION    TMP1(irun),   TMP2(irun)  
       DIMENSION    TMP3(irun),  ITMP1(irun)  
       DIMENSION   ITMP2(irun)  
       dimension qqcolmin(irun),qqcolmax(irun),levpbl(irun)  
   
 C Dynamic Variables  
 C -----------------  
       DIMENSION ADZ1   (irun,nlev  ), DZ1TMP(irun,nlev)  
       DIMENSION DZ3    (irun,nlev  ), TEMP  (irun,nlev)  
       DIMENSION DV     (irun,nlev  ), DTHV  (irun,nlev)  
       DIMENSION DPK    (irun,nlev  ), STRT  (irun,nlev)  
       DIMENSION DW2    (irun,nlev  ), RI    (irun,nlev)  
       DIMENSION RHOZPK (irun,nlev  ), Q     (irun,nlev)  
       DIMENSION RIINIT (irun,nlev  ), DU    (irun,nlev)  
       DIMENSION QQINIT (irun,nlev  ), RHOKDZ(irun,nlev)  
       DIMENSION RHODZ2 (irun,nlev  )  
       REAL      KM     (irun,nlev  ),     KH(irun,nlev)  
   
 C Dynamic Variables  
 C -----------------  
       DIMENSION DELTH  (irun,nlev+1), DELSH (irun,nlev+1)  
       DIMENSION FLXFAC (irun,nlev+1), DTHG  (irun,nlev+1)  
       DIMENSION FLXFPK (irun,nlev+1)  
   
 C Dynamic Variables  
 C -----------------  
       DIMENSION ADZ2   (irun,nlev-1), RHODZ1(irun,nlev-1)  
       DIMENSION VKZE   (irun,nlev-1), VKZM  (irun,nlev-1)  
       DIMENSION XL     (irun,nlev-1), QXLM  (irun,nlev-1)  
       DIMENSION QQE    (irun,nlev-1), QE    (irun,nlev-1)  
       DIMENSION P3     (irun,nlev-1), XQ    (irun,nlev-1)  
       DIMENSION XLDIAG (irun,nlev-1), FLXFCE(irun,nlev-1)  
1833    
1834        LOGICAL FIRST,LAST        LOGICAL FIRST,LAST
1835        DIMENSION IBITSTB(irun,nlev),IBITSTR(irun),INTQ(irun,nlev)        integer IBITSTB(irun,nlev),INTQ(irun,nlev)
1836    
1837  C arrays for use by moist bouyancy calculation  C arrays for use by moist bouyancy calculation
1838  C -----------------  C -----------------
1839        real TL(irun,NLEV),DTH(irun,NLEV)        _RL TL(irun,NLEV),DTH(irun,NLEV)
1840        real DSH(irun,NLEV)        _RL DSH(irun,NLEV)
1841        real SHL(irun,NLEV)        _RL SHL(irun,NLEV)
1842        real AA(irun,NLEV),BB(irun,NLEV),SSDEV(irun,NLEV)        _RL AA(irun,NLEV),BB(irun,NLEV),SSDEV(irun,NLEV)
1843        real ARG(irun,NLEV),XXZETA(irun),QBYU(irun)        _RL ARG(irun,NLEV),XXZETA(irun),QBYU(irun)
1844        real SVAR(irun,NLEV),Q1M(irun,NLEV)        _RL SVAR(irun,NLEV),Q1M(irun,NLEV)
1845        real FCC(irun,NLEV)        _RL FCC(irun,NLEV)
1846        real BETAT(irun,NLEV),BETAW(irun,NLEV)        _RL BETAT(irun,NLEV),BETAW(irun,NLEV)
1847        real BETAL(irun,NLEV),BETAT1(irun,NLEV)        _RL BETAL(irun,NLEV),BETAT1(irun,NLEV)
1848        real BETAW1(irun,NLEV),SBAR(irun,NLEV)        _RL BETAW1(irun,NLEV),SBAR(irun,NLEV)
1849        real SHSAT(irun,NLEV)        _RL SHSAT(irun,NLEV)
       real TEMPOR(irun,NLEV)  
1850    
1851  C Some space for variables to be used in called routines  C Some space for variables to be used in called routines
1852        logical LWATER        logical LWATER
1853        integer IVBITRIB(irun)        integer IVBITRIB(irun)
1854        DIMENSION VHZ(irun)        _RL VHZ(irun)
1855        DIMENSION VH0(irun)        _RL VH0(irun)
1856        DIMENSION VPSIM(irun),VAPSIM(irun)        _RL VPSIM(irun),VAPSIM(irun)
1857        DIMENSION VPSIG(irun),VPSIHG(irun)        _RL VPSIG(irun),VPSIHG(irun)
1858        DIMENSION VTEMP(irun),VDZETA(irun)        _RL VTEMP(irun),VDZETA(irun)
1859        DIMENSION VDZ0(irun),VDPSIM(irun)        _RL VDZ0(irun),VDPSIM(irun)
1860        DIMENSION VDPSIH(irun),VZH(irun)        _RL VDPSIH(irun),VZH(irun)
1861        DIMENSION VXX0(irun),VYY0(irun)        _RL VXX0(irun),VYY0(irun)
1862        DIMENSION VAPSIHG(irun),VRIB1(irun),VWS1(irun)        _RL VAPSIHG(irun),VRIB1(irun),VWS1(irun)
1863        DIMENSION VPSIH(irun),VZETAL(irun)        _RL VPSIH(irun),VZETAL(irun)
1864        DIMENSION VZ0L(irun),VPSIH2(irun)        _RL VZ0L(irun),VPSIH2(irun)
1865        DIMENSION VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)        _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
1866        DIMENSION VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)        _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
1867        DIMENSION VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)        _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
1868        DIMENSION VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)        _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
1869        DIMENSION VXNUM(irun),VDZETA1(irun),VDZETA2(irun)        _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
1870        DIMENSION VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)        _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
1871        DIMENSION VDPSIMC(irun),VDPSIHC(irun)        _RL VDPSIMC(irun),VDPSIHC(irun)
1872        integer types(irun)  
1873        character*40 name        _RL DZITRP(irun,nlev-1),STBFCN(irun,nlev)
1874          _RL XL0(irun,nlev),Q1(irun,nlev-1)
1875        DIMENSION DZITRP(irun,nlev-1), STBFCN(irun,nlev)        _RL WRKIT1(irun,nlev-1)
1876        DIMENSION    XL0(irun,nlev),       Q1(irun,nlev-1)        _RL WRKIT2(irun,nlev-1)
1877        DIMENSION WRKIT1(irun,nlev-1)        _RL WRKIT3(irun,nlev-1)
1878        DIMENSION WRKIT2(irun,nlev-1)        _RL WRKIT4(irun,nlev-1)
       DIMENSION WRKIT3(irun,nlev-1)  
       DIMENSION WRKIT4(irun,nlev-1)  
1879        INTEGER INT1(irun,nlev), INT2(irun,nlev-1)        INTEGER INT1(irun,nlev), INT2(irun,nlev-1)
1880    
1881        real vrt1con,pi,rsq2pi,p5sr,clh        _RL vrt1con,pi,rsq2pi,p5sr,clh,vk,rvk,aitr,gbycp,fac1,fac2
1882        integer nt        _RL getcon,dum,errf
1883          integer istnlv,nlevm1,nlevm2,nlevml,nlevp1,istnm1,istnm2,istnp1
1884          integer istnml,istnmq,istlmq,nlevmq
1885          integer i,iter,init,n,nt,LL,L,Lp,Lp1,lmin,lminq,lminq1,ibit
1886    
1887        vk = getcon('VON KARMAN')        vk = getcon('VON KARMAN')
1888        rvk = 1./vk        rvk = 1./vk
# Line 1722  C Some space for variables to be used in Line 1902  C Some space for variables to be used in
1902        P5SR = 0.5**0.5        P5SR = 0.5**0.5
1903        CLH = GETCON('LATENT HEAT COND') / CP        CLH = GETCON('LATENT HEAT COND') / CP
1904    
 C CALL POINT BY POINT DIAGNOSTIC ROUTINE FOR 'BEFORE' VALUES  
 C ----------------------------------------------------------  
       IF(TPROF) THEN  
       do i = 1,irun  
        types(i) = 1  
       enddo  
       name = 'mid pressure'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,pl)  
       name = 'edge pressure'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev+1,name,ple)  
       name = 'mid p ** kappa'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,plke)  
       name = 'edge p ** kappa'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev+1,name,plke)  
       name = 'theta before'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev+1,name,th)  
       name = 'theta virtual before'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev+1,name,thv)  
       name = 'q before'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev+1,name,sh)  
       name = 'u wind before'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,u)  
       name = 'v wind before'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,V)  
       name = 'heat cap ground'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,hcapg)  
       name = 'latent heat at ground'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,clhg)  
       name = 'net surface rad'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,radflx)  
       name = 'background heat transfer'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,khbg)  
       ENDIF  
1905  C SET INITIAL NUMBER OF ITERATIONS OF SFCFLX  C SET INITIAL NUMBER OF ITERATIONS OF SFCFLX
1906  C ------------------------------------------  C ------------------------------------------
1907        N = 6        N = 6
# Line 2314  C ---------------------------------- Line 2461  C ----------------------------------
2461       1      * QXLM(I,LMINQ1)       1      * QXLM(I,LMINQ1)
2462  9306  CONTINUE  9306  CONTINUE
2463         CALL TRBDIF(QQ(1,LMINQ1),P3(1,LMINQ1),RHOKDZ(1,LMINQ1),         CALL TRBDIF(QQ(1,LMINQ1),P3(1,LMINQ1),RHOKDZ(1,LMINQ1),
2464       1       FLXFCE(1,LMINQ1),DTHS,DELTHS,NLEVMQ,1,1.0E-20,irun)       1       FLXFCE(1,LMINQ1),DTHS,DELTHS,NLEVMQ,1,1.0 _d -20,irun)
2465  C  C
2466  C SAVE KH BEFORE ADDING DIMENSIONS FOR USE BY MOIST BOUYANCY CALCULATION  C SAVE KH BEFORE ADDING DIMENSIONS FOR USE BY MOIST BOUYANCY CALCULATION
2467  C  C
# Line 2355  C Line 2502  C
2502        DO 9132 I =1,irun        DO 9132 I =1,irun
2503        DELTH(I,NLEVP1) = 1.        DELTH(I,NLEVP1) = 1.
2504  9132  CONTINUE  9132  CONTINUE
2505         CALL TRBDIF(TH,DELTH,TEMP,FLXFPK,DTHS,DELTHS,NLEV,2,0.,irun)         CALL TRBDIF(TH,DELTH,TEMP,FLXFPK,DTHS,DELTHS,NLEV,2,0. _d 0,irun)
2506        do i = 1,irun        do i = 1,irun
2507         hsturb(i) = -1.* dths(i)         hsturb(i) = -1.* dths(i)
2508         dhsdtc(i) = -1.* delths(i)         dhsdtc(i) = -1.* delths(i)
# Line 2381  C Line 2528  C
2528        DELSH(I,NLEVP1) = 1.        DELSH(I,NLEVP1) = 1.
2529  9140  CONTINUE  9140  CONTINUE
2530    
2531         CALL TRBDIF(SH,DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,NLEV,2,0.,irun)         CALL TRBDIF(SH,DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,NLEV,
2532         .                                     2,0. _d 0,irun)
2533        do i = 1,irun        do i = 1,irun
2534         eturb(i) = -1.* dthl(i)         eturb(i) = -1.* dthl(i)
2535         dedqa(i) = -1.* delthl(i)         dedqa(i) = -1.* delthl(i)
# Line 2404  C Line 2552  C
2552        rhokdz(i,nlev) = 0.0              rhokdz(i,nlev) = 0.0      
2553        enddo        enddo
2554    
2555        do nt = 1,ntrace  c     do nt = 1,ntrace
2556        do  i = 1,irun  c     do  i = 1,irun
2557        tracers(i,nlev+1,nt) = tracers(i,nlev,nt)  c     tracers(i,nlev+1,nt) = tracers(i,nlev,nt)
2558        enddo  c     enddo
2559        CALL TRBDIF(tracers(1,1,nt),DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,  c     CALL TRBDIF(tracers(1,1,nt),DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,
2560       .                                                  NLEV,4,0.,irun)  c    .                                         NLEV,4,0. _d 0,irun)
2561        enddo  c     enddo
2562  C  C
2563  C   CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: U AND V  C   CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: U AND V
2564  C  C
2565        DO 9172 I =1,ISTNLV        DO 9172 I =1,ISTNLV
2566        RHOKDZ(I,1) = RHODZ2(I,1) * KM(I,1)        RHOKDZ(I,1) = RHODZ2(I,1) * KM(I,1)
2567  9172  CONTINUE  9172  CONTINUE
2568         CALL TRBDIF(U,V,RHOKDZ,FLXFAC,DTHS,DELTHS,NLEV,3,0.,irun)         CALL TRBDIF(U,V,RHOKDZ,FLXFAC,DTHS,DELTHS,NLEV,3,0. _d 0,irun)
2569  C ( FILL DIAGNOSTIC ARRAYS IF REQUIRED )  C ( FILL DIAGNOSTIC ARRAYS IF REQUIRED )
2570        DO 9174 I =1,ISTNLV        DO 9174 I =1,ISTNLV
2571        WU(I,1) = WU(I,1) + RHOKDZ(I,1) * ( U(I,2) - U(I,1) )        WU(I,1) = WU(I,1) + RHOKDZ(I,1) * ( U(I,2) - U(I,1) )
# Line 2514  C Line 2662  C
2662         WV(I,1) =  WV(I,1) * AITR         WV(I,1) =  WV(I,1) * AITR
2663  9194  CONTINUE  9194  CONTINUE
2664  C  C
 C  IF TPROF, CALL POINT BY POINT DIAGNOSTIC ROUTINE FOR 'AFTER' VALUES  
 C  
        IF(TPROF)THEN  
       name = 'tke before'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,qqinit)  
       name = 'richardson number before'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,riinit)  
       name = 'theta after'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev+1,name,th)  
       name = 'theta virtual after'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev+1,name,thv)  
       name = 'q after'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev+1,name,sh)  
       name = 'u wind after'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,u)  
       name = 'v wind after'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,V)  
       name = 'tke after'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,qq)  
       name = 'richardson number after'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,ri)  
       name = 'trb u flux'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,wu)  
       name = 'trb v flux'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,wv)  
       name = 'trb t flux'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,wt)  
       name = 'trb q flux'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,wsh)  
       name = 'eddy coef mom'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,eu)  
       name = 'eddy coef heat'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,et)  
       name = 'length scale'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,xldiag)  
       name = 'layer heights'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,nlev,name,dz1tmp)  
       name = 'q ground'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,shg)  
       name = 'rib initial'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,ribin)  
       name = 'cu initial'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,cuin)  
       name = 'ct initial'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,ctin)  
       name = 'ustar initial'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,ustarin)  
       name = 'rho sfc initial'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,rhosin)  
       name = 'z0 initial'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,z0in)  
       name = 'zeta initial'  
       call pntprf(irun,irun,lons,lats,types,nymd,nhms,1,name,zetain)  
        ENDIF  
2665        RETURN        RETURN
2666        END        END
2667         SUBROUTINE SFCFLX(NN,VUS,VVS,VTHV1,VTHV2,VTH1,VTH2,VSH1,         SUBROUTINE SFCFLX(NN,VUS,VVS,VTHV1,VTHV2,VTH1,VTH2,VSH1,
# Line 2622  C    CU            -         MOMENTUM TR Line 2716  C    CU            -         MOMENTUM TR
2716  C    CT            -         HEAT TRANSPORT COEFFICIENT  C    CT            -         HEAT TRANSPORT COEFFICIENT
2717  C  C
2718  C**********************************************************************  C**********************************************************************
2719  C        implicit none
2720    
2721    C Argument List Declarations
2722          integer nn,n,irun
2723          _RL aitr,cp,rgas,undef
2724          _RL VUS(IRUN),VVS(IRUN),VTHV1(IRUN),VTHV2(IRUN)
2725          _RL VTH1(IRUN),VTH2(IRUN),VSH1(IRUN),VSH2(IRUN)
2726          _RL VPK(IRUN),VPKE(IRUN),VPE(IRUN)
2727          _RL VZ0(IRUN),VHS(IRUN),VAHS(IRUN)
2728          integer IVWATER(IRUN)
2729          LOGICAL FIRST,LAST
2730          _RL VRHO(IRUN),VRHOZPK(IRUN)
2731          _RL VKM(IRUN),VKH(IRUN),VUSTAR(IRUN),VXX(IRUN)
2732          _RL VYY(IRUN),VCU(IRUN),VCT(IRUN),VRIB(IRUN)
2733          _RL VZETA(IRUN),VWS(IRUN)
2734          _RL stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)
2735          _RL stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)
2736          LOGICAL LWATER
2737          integer IVBITRIB(irun)
2738          _RL VHZ(irun),VPSIM(irun),VAPSIM(irun),VPSIG(irun),VPSIHG(irun)
2739          _RL VTEMP(irun),VDZETA(irun),VDZ0(irun),VDPSIM(irun)
2740          _RL VDPSIH(irun),VZH(irun),VXX0(irun),VYY0(irun)
2741          _RL VAPSIHG(irun),VRIB1(irun),VWS1(irun)
2742          _RL VPSIH(irun),VZETAL(irun),VZ0L(irun),VPSIH2(irun),VH0(irun)
2743          _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
2744          _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
2745          _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
2746          _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
2747          _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
2748          _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
2749          _RL VDPSIMC(irun),VDPSIHC(irun)
2750    
2751    C Local Variables
2752          _RL USTMX3,USTZ0S,Z0MIN,H0BYZ0,USTH0S,H0VEG,Z0VEGM,PRFAC
2753          _RL XPFAC,DIFSQT
2754        PARAMETER ( USTMX3 =   0.0632456)        PARAMETER ( USTMX3 =   0.0632456)
2755        PARAMETER ( USTZ0S =   0.2030325E-5)        PARAMETER ( USTZ0S =   0.2030325E-5)
2756        PARAMETER ( Z0MIN  =  USTZ0S/USTMX3)        PARAMETER ( Z0MIN  =  USTZ0S/USTMX3)
# Line 2635  C Line 2762  C
2762        PARAMETER ( XPFAC  = .55        )          PARAMETER ( XPFAC  = .55        )  
2763        PARAMETER ( DIFSQT  = 3.872983E-3)        PARAMETER ( DIFSQT  = 3.872983E-3)
2764    
2765         DIMENSION VUS(IRUN),VVS(IRUN),VTHV1(IRUN),VTHV2(IRUN)        _RL psihdiag(irun),psimdiag(irun)
2766         DIMENSION VTH1(IRUN),VTH2(IRUN),VSH1(IRUN),VSH2(IRUN)        _RL getcon,vk,rvk,vk2,bmdl
2767         DIMENSION VPK(IRUN),VPKE(IRUN),VPE(IRUN)        integer iwater,itype
2768         DIMENSION VZ0(IRUN),IVWATER(IRUN),VHS(IRUN),VAHS(IRUN)        integer i,iter
        DIMENSION VRHO(IRUN),VRHOZPK(IRUN)  
        DIMENSION VKM(IRUN),VKH(IRUN),VUSTAR(IRUN),VXX(IRUN)  
        DIMENSION VYY(IRUN),VCU(IRUN),VCT(IRUN),VRIB(IRUN)  
        DIMENSION VZETA(IRUN),VWS(IRUN)  
        dimension stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)  
        dimension stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)  
        LOGICAL FIRST,LAST  
        LOGICAL LWATER  
        integer IVBITRIB(irun)  
 C  
        DIMENSION VHZ(irun)  
        DIMENSION VH0(irun)  
        DIMENSION VPSIM(irun),VAPSIM(irun)  
        DIMENSION VPSIG(irun),VPSIHG(irun)  
        DIMENSION VTEMP(irun),VDZETA(irun)  
        DIMENSION VDZ0(irun),VDPSIM(irun)  
        DIMENSION VDPSIH(irun),VZH(irun)  
        DIMENSION VXX0(irun),VYY0(irun)  
        DIMENSION VAPSIHG(irun),VRIB1(irun),VWS1(irun)  
        DIMENSION VPSIH(irun),VZETAL(irun)  
        DIMENSION VZ0L(irun),VPSIH2(irun)  
        DIMENSION VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)  
        DIMENSION VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)  
        DIMENSION VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)  
        DIMENSION VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)  
        DIMENSION VXNUM(irun),VDZETA1(irun),VDZETA2(irun)  
        DIMENSION VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)  
        DIMENSION VDPSIMC(irun),VDPSIHC(irun)  
   
        dimension psihdiag(irun),psimdiag(irun)  
2769  C  C
2770        vk = getcon('VON KARMAN')        vk = getcon('VON KARMAN')
2771        rvk = 1./vk        rvk = 1./vk
# Line 3005  C                  (IFLAG=3), OR BOTH (I Line 3102  C                  (IFLAG=3), OR BOTH (I
3102  C     N     -  LENGTH OF VECTOR TO BE SOLVED  C     N     -  LENGTH OF VECTOR TO BE SOLVED
3103  C  C
3104  C**********************************************************************  C**********************************************************************
3105  C        implicit none
3106        DIMENSION PHIM(N),PHIH(N),Z(N)  
3107        DIMENSION INT72(N), INTMAX(N)  C Argument List Declarations
3108        DIMENSION ZSTAR(N),I1(N),I2(N)        integer n,iflag
3109        DIMENSION E1(N),E2(N),TEMP1(N)        _RL PHIM(N),PHIH(N),Z(N)
3110  C  
3111        DIMENSION PHIM0(385),ZLINM1(75),ZLINM2(75),ZLINM3(36)  C Local Variables
3112        DIMENSION ZLOGM1(74),ZLOGM2(75),ZLOGM3(50)        integer I1(N),I2(N)
3113        DIMENSION PHIH0(385),ZLINH1(75),ZLINH2(75),ZLINH3(36)        _RL ZSTAR(N),E1(N),E2(N),TEMP1(N)
3114        DIMENSION ZLOGH1(74),ZLOGH2(75),ZLOGH3(50)  C
3115          _RL PHIM0(385),ZLINM1(75),ZLINM2(75),ZLINM3(36)
3116          _RL ZLOGM1(74),ZLOGM2(75),ZLOGM3(50)
3117          _RL PHIH0(385),ZLINH1(75),ZLINH2(75),ZLINH3(36)
3118          _RL ZLOGH1(74),ZLOGH2(75),ZLOGH3(50)
3119        EQUIVALENCE (PHIM0(1),ZLINM1(1)),(PHIM0(76),ZLINM2(1))        EQUIVALENCE (PHIM0(1),ZLINM1(1)),(PHIM0(76),ZLINM2(1))
3120        EQUIVALENCE (PHIM0(151),ZLINM3(1))        EQUIVALENCE (PHIM0(151),ZLINM3(1))
3121        EQUIVALENCE (PHIM0(187),ZLOGM1(1)),(PHIM0(261),ZLOGM2(1))        EQUIVALENCE (PHIM0(187),ZLOGM1(1)),(PHIM0(261),ZLOGM2(1))
# Line 3193  C Line 3294  C
3294       &  0.664746,0.663985,0.663227,0.662473,0.661723,       &  0.664746,0.663985,0.663227,0.662473,0.661723,
3295       &  0.660977,0.660234,0.659495,0.658759,0.658027,       &  0.660977,0.660234,0.659495,0.658759,0.658027,
3296       &  0.657298/       &  0.657298/
3297    
3298          integer ibit1,ibit2,i
3299  C  C
3300        IBIT1 = 0        IBIT1 = 0
3301        IBIT2 = 0        IBIT2 = 0
# Line 3299  C  SUBPROGRAMS NEEDED Line 3402  C  SUBPROGRAMS NEEDED
3402  C     PHI  -  COMPUTES SIMILARITY FUNCTION FOR MOMENTUM AND SCALARS  C     PHI  -  COMPUTES SIMILARITY FUNCTION FOR MOMENTUM AND SCALARS
3403  C  C
3404  C**********************************************************************  C**********************************************************************
3405  C        implicit none
3406  C  
3407    C Argument List Declarations
3408          integer irun,iflag
3409          _RL VZZ(IRUN),VZH(IRUN),VPSIM(IRUN),VPSIH(IRUN),
3410         1     VX(IRUN),VXS(IRUN),VY(IRUN),VYS(IRUN)
3411    
3412    C Local Variables
3413          _RL ZWM,RZWM,Z0M,ZCM,RZCM,CM1,CM2,CM6,CM7,CM8ARG,YCM
3414        PARAMETER ( ZWM     =    1.    )        PARAMETER ( ZWM     =    1.    )
3415        PARAMETER ( RZWM    =  1./ZWM  )        PARAMETER ( RZWM    =  1./ZWM  )
3416        PARAMETER ( Z0M     =    0.2    )        PARAMETER ( Z0M     =    0.2    )
# Line 3313  C Line 3423  C
3423        PARAMETER ( CM8ARG  =  CM7*ZCM*RZWM / (CM2+ZCM)  )        PARAMETER ( CM8ARG  =  CM7*ZCM*RZWM / (CM2+ZCM)  )
3424        PARAMETER ( YCM     =  6. / ( 1. + 6.*CM1*ZCM )  )        PARAMETER ( YCM     =  6. / ( 1. + 6.*CM1*ZCM )  )
3425    
3426        DIMENSION VZZ(IRUN),VZH(IRUN),VPSIM(IRUN),VPSIH(IRUN),        integer INTSTB(irun),INTZ0(irun)
3427       1     VX(IRUN),VXS(IRUN),VY(IRUN),VYS(IRUN)        _RL ZZ0(irun),Z(irun),Z2(irun),Z1(irun),Z0(irun)
3428          _RL X0(irun),X1(irun),Y0(irun),Y1(irun)
3429        DIMENSION INTSTB(irun),INT72(irun),INTZ0(irun)        _RL PSI2(irun),TEMP(irun)
3430        DIMENSION  ZZ0(irun),Z(irun),Z2(irun),Z1(irun),Z0(irun)        _RL HZ(irun),ARG0(irun),ARG1(irun),DX(irun)
3431        DIMENSION  X0(irun),X1(irun),Y0(irun),Y1(irun)        _RL X0NUM(irun),X1NUM(irun),X0DEN(irun)
3432        DIMENSION  XX(irun),XXS(irun),YY(irun),YYS(irun)        _RL X1DEN(irun),Y1DEN(irun),Z2ZWM(irun)
3433        DIMENSION  PSIMM(irun),PSIHH(irun)        _RL cm3,cm4,cm5,cm8
3434        DIMENSION  PSI2(irun),TEMP(irun)        integer ibit,indx
3435        DIMENSION  HZ(irun),ARG0(irun),ARG1(irun),DX(irun)        integer i
       DIMENSION  X0NUM(irun),X1NUM(irun),X0DEN(irun)  
       DIMENSION  X1DEN(irun),Y1DEN(irun),Z2ZWM(irun)  
3436  C  C
3437        CM3 =   sqrt( 0.2/CM1-0.01 )        CM3 =   sqrt( 0.2/CM1-0.01 )
3438        CM4 =   1./CM3        CM4 =   1./CM3
# Line 3358  C ************************************** Line 3466  C **************************************
3466  C  C
3467        IF(IBIT.LE.0)  GO TO 100        IF(IBIT.LE.0)  GO TO 100
3468  C  C
3469        INDEX = 0        indx = 0
3470        DO 9002 I = 1,IRUN        DO 9002 I = 1,IRUN
3471         IF (INTSTB(I).EQ.1)THEN         IF (INTSTB(I).EQ.1)THEN
3472          INDEX = INDEX + 1          indx = indx + 1
3473          Z(INDEX) = VZZ(I)          Z(indx) = VZZ(I)
3474          Z0(INDEX) = ZZ0(I)          Z0(indx) = ZZ0(I)
3475         ENDIF         ENDIF
3476   9002 CONTINUE   9002 CONTINUE
3477  C  C
# Line 3401  C Line 3509  C
3509         PSI2(I) = PSI2(I) + DX(I)         PSI2(I) = PSI2(I) + DX(I)
3510   9006 CONTINUE   9006 CONTINUE
3511  C  C
3512        INDEX = 0        indx = 0
3513        DO 9008 I = 1,IRUN        DO 9008 I = 1,IRUN
3514         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3515          INDEX = INDEX + 1          indx = indx + 1
3516          VPSIM(I) = PSI2(INDEX)          VPSIM(I) = PSI2(indx)
3517          VX(I) = X1(INDEX)          VX(I) = X1(indx)
3518          VXS(I) = X0(INDEX)          VXS(I) = X0(indx)
3519         ENDIF         ENDIF
3520   9008 CONTINUE   9008 CONTINUE
3521  C  C
# Line 3434  C Line 3542  C
3542         PSI2(I) = PSI2(I) - Y1(I) + Y0(I)         PSI2(I) = PSI2(I) - Y1(I) + Y0(I)
3543   9010 CONTINUE   9010 CONTINUE
3544  C  C
3545        INDEX = 0        indx = 0
3546        DO 9012 I = 1,IRUN        DO 9012 I = 1,IRUN
3547         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3548         INDEX = INDEX + 1         indx = indx + 1
3549         VPSIH(I) = PSI2(INDEX)         VPSIH(I) = PSI2(indx)
3550         VY(I) = Y1(INDEX)         VY(I) = Y1(indx)
3551         VYS(I) = Y0(INDEX)         VYS(I) = Y0(indx)
3552         ENDIF         ENDIF
3553   9012 CONTINUE   9012 CONTINUE
3554  C  C
# Line 3463  C Line 3571  C
3571         ENDIF         ENDIF
3572   9014 CONTINUE   9014 CONTINUE
3573        IF(IBIT.LE.0)  GO TO 300        IF(IBIT.LE.0)  GO TO 300
3574        INDEX = 0        indx = 0
3575  #if CRAY  #ifdef CRAY
3576  CDIR$ NOVECTOR  CDIR$ NOVECTOR
3577  #endif  #endif
3578        DO 9016 I = 1,IRUN        DO 9016 I = 1,IRUN
3579         IF (INTSTB(I).EQ.1)THEN         IF (INTSTB(I).EQ.1)THEN
3580          INDEX = INDEX + 1          indx = indx + 1
3581          Z(INDEX) = VZZ(I)          Z(indx) = VZZ(I)
3582          Z0(INDEX) = ZZ0(I)          Z0(indx) = ZZ0(I)
3583          ARG1(INDEX) = VZH(I)          ARG1(indx) = VZH(I)
3584         ENDIF         ENDIF
3585   9016 CONTINUE   9016 CONTINUE
3586  #if CRAY  #ifdef CRAY
3587  CDIR$ VECTOR  CDIR$ VECTOR
3588  #endif  #endif
3589    
# Line 3529  C Line 3637  C
3637         PSI2(I) = TEMP(I) + CM6 * ARG1(I)         PSI2(I) = TEMP(I) + CM6 * ARG1(I)
3638   9020 CONTINUE   9020 CONTINUE
3639  C  C
3640        INDEX = 0        indx = 0
3641        DO 9030 I = 1,IRUN        DO 9030 I = 1,IRUN
3642         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3643         INDEX = INDEX + 1         indx = indx + 1
3644         VPSIM(I) = PSI2(INDEX)         VPSIM(I) = PSI2(indx)
3645         VX(I) = X1(INDEX)         VX(I) = X1(indx)
3646         VXS(I) = X0(INDEX)         VXS(I) = X0(indx)
3647         ENDIF         ENDIF
3648   9030 CONTINUE   9030 CONTINUE
3649  C  C
# Line 3561  C Line 3669  C
3669         PSI2(I) = TEMP(I) + ARG0(I) * ARG1(I)         PSI2(I) = TEMP(I) + ARG0(I) * ARG1(I)
3670   9024 CONTINUE   9024 CONTINUE
3671  C  C
3672        INDEX = 0        indx = 0
3673        DO 9026 I = 1,IRUN        DO 9026 I = 1,IRUN
3674         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3675         INDEX = INDEX + 1         indx = indx + 1
3676         VPSIH(I) = PSI2(INDEX)         VPSIH(I) = PSI2(indx)
3677         VY(I) = Y1(INDEX)         VY(I) = Y1(indx)
3678         VYS(I) = X0(INDEX)         VYS(I) = X0(indx)
3679         ENDIF         ENDIF
3680   9026 CONTINUE   9026 CONTINUE
3681  C  C
# Line 3613  C Line 3721  C
3721  C     TRBITP -  INTERPOLATES TO HEIGHT WHERE RI = RICR  C     TRBITP -  INTERPOLATES TO HEIGHT WHERE RI = RICR
3722  C  C
3723  C**********************************************************************  C**********************************************************************
3724  C        implicit none
3725  C  
3726    C Argument List Declarations
3727          integer irun,nlev,init,lmin,lminq,lminq1
3728          _RL cp
3729          _RL STRT(irun,NLEV),DW2(irun,NLEV),DZ3(irun,NLEV)
3730          _RL Q(irun,NLEV),VKZM(irun,NLEV-1),VKZE(irun,NLEV-1)
3731          _RL DTHV(irun,NLEV),DPK(irun,NLEV),DU(irun,NLEV)
3732          _RL DV(irun,NLEV)
3733          _RL QXLM(irun,NLEV-1),XL(irun,NLEV-1)
3734          _RL DZITRP(irun,nlev-1),STBFCN(irun,nlev)
3735          _RL XL0(irun,nlev),Q1(irun,nlev-1)
3736          _RL WRKIT1(irun,nlev-1),WRKIT2(irun,nlev-1)
3737          _RL WRKIT3(irun,nlev-1)
3738          _RL WRKIT4(irun,nlev-1)
3739          INTEGER INT1(irun,nlev), INT2(irun,nlev-1)
3740    
3741    C Local Variables
3742          _RL rf1,rf2,e5,d4,d1,rfc,ricr,alpha,dzcnv,xl0cnv,xl0min
3743          _RL clmt,clmt53
3744        PARAMETER ( RF1     = 0.2340678 )        PARAMETER ( RF1     = 0.2340678 )
3745        PARAMETER ( RF2     = 0.2231172 )        PARAMETER ( RF2     = 0.2231172 )
3746        PARAMETER ( E5      = 49.66     )          PARAMETER ( E5      = 49.66     )  
# Line 3628  C Line 3754  C
3754        PARAMETER ( XL0MIN  = 1.        )        PARAMETER ( XL0MIN  = 1.        )
3755        PARAMETER ( CLMT    = 0.23      )          PARAMETER ( CLMT    = 0.23      )  
3756        PARAMETER ( CLMT53  = 5. * CLMT / 3. )        PARAMETER ( CLMT53  = 5. * CLMT / 3. )
3757    
3758          integer ibit,nlevm1,nlevp1,istnlv,istnm1,nlevml,istnml,Lp1
3759        DIMENSION STRT(irun,NLEV),   DW2(irun,NLEV), DZ3(irun,NLEV)        integer istnmq,istlmq,lminp,lm1,lmin1
3760        DIMENSION DTHV(irun,NLEV),   DPK(irun,NLEV),  DU(irun,NLEV)        integer i,L,LL
       DIMENSION   DV(irun,NLEV),     Q(irun,NLEV)  
       DIMENSION VKZM(irun,NLEV-1),VKZE(irun,NLEV-1)  
       DIMENSION QXLM(irun,NLEV-1),  XL(irun,NLEV-1)  
       DIMENSION DZITRP(irun,nlev-1), STBFCN(irun,nlev)  
       DIMENSION    XL0(irun,nlev),       Q1(irun,nlev-1)  
       DIMENSION WRKIT1(irun,nlev-1)  
       DIMENSION WRKIT2(irun,nlev-1)  
       DIMENSION WRKIT3(irun,nlev-1)  
       DIMENSION WRKIT4(irun,nlev-1)  
 C  
       INTEGER INT1(irun,nlev), INT2(irun,nlev-1)  
3761  C  C
3762        NLEVM1 = NLEV - 1        NLEVM1 = NLEV - 1
3763        NLEVP1 = NLEV + 1        NLEVP1 = NLEV + 1
# Line 3846  C     ------- Line 3961  C     -------
3961  C    DZITRP        -         INTERPOLATION COEFFICIENT  C    DZITRP        -         INTERPOLATION COEFFICIENT
3962  C  C
3963  C**********************************************************************  C**********************************************************************
3964  C        implicit none
3965  C  
3966    C Argument List Declarations
3967          integer irun,nlev
3968          _RL cp
3969          _RL STBFCN(irun,NLEV+1)
3970          integer INTCHG(irun,NLEV)
3971          _RL DTHV(irun,NLEV+1),DPK(irun,NLEV+1)
3972          _RL DU(irun,NLEV+1),DV(irun,NLEV+1)
3973          _RL DZITRP(irun,NLEV+1)
3974          _RL AAA(irun,NLEV),BBB(irun,NLEV)
3975          _RL CCC(irun,NLEV),DDD(irun,NLEV)
3976    
3977    C Local Variables
3978          _RL rf1,rf2,e5,d4,d1,rfc,ricr
3979        PARAMETER ( RF1     = 0.2340678 )        PARAMETER ( RF1     = 0.2340678 )
3980        PARAMETER ( RF2     = 0.2231172 )        PARAMETER ( RF2     = 0.2231172 )
3981        PARAMETER ( E5      = 49.66     )          PARAMETER ( E5      = 49.66     )  
# Line 3856  C Line 3984  C
3984        PARAMETER ( RFC     = 0.1912323 )        PARAMETER ( RFC     = 0.1912323 )
3985        PARAMETER ( RICR    = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 )  )        PARAMETER ( RICR    = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 )  )
3986    
3987        DIMENSION STBFCN(irun,NLEV+1), INTCHG(irun,NLEV)        integer istnlv
3988        DIMENSION DTHV  (irun,NLEV+1), DPK   (irun,NLEV+1)        integer i
       DIMENSION DU    (irun,NLEV+1), DV    (irun,NLEV+1)  
       DIMENSION DZITRP(irun,NLEV+1)  
       DIMENSION AAA   (irun,NLEV), BBB   (irun,NLEV)  
       DIMENSION CCC   (irun,NLEV), DDD   (irun,NLEV)  
3989  C  C
3990  C *********************************************************************  C *********************************************************************
3991  C ****           QUADRATIC INTERPOLATION OF RI TO RICR VIA          ***  C ****           QUADRATIC INTERPOLATION OF RI TO RICR VIA          ***
# Line 3924  C    QQE           -         EQUILIBRIUM Line 4048  C    QQE           -         EQUILIBRIUM
4048  C    BITSTB        -         BIT '1' WHERE QE GREATER THAN ZERO  C    BITSTB        -         BIT '1' WHERE QE GREATER THAN ZERO
4049  C  C
4050  C**********************************************************************  C**********************************************************************
4051  C        implicit none
4052  C  
4053    C Argument List Declarations
4054          integer nlev,nlay,irun
4055          _RL RI(irun,NLEV),STRT(irun,NLEV),DW2(irun,NLEV)
4056          _RL XL(irun,NLEV),ZKM(irun,NLEV),ZKH(irun,NLEV)
4057          _RL QE(irun,NLEV),QQE(irun,NLEV)
4058          INTEGER INTSTB(irun,nlev)
4059          _RL EE(irun,nlay-1),RF(irun,nlay-1)
4060    
4061    C Local Variables
4062          _RL b1,b2,d3,rf1,rf2,d3b2,d2,e5,d4,d1,d1half,d2half
4063          _RL rfc,ricr,ch,cm
4064        PARAMETER ( B1      =   16.6    )          PARAMETER ( B1      =   16.6    )  
4065        PARAMETER ( B2      =   10.1    )        PARAMETER ( B2      =   10.1    )
4066        PARAMETER ( D3      = 0.29397643 )        PARAMETER ( D3      = 0.29397643 )
# Line 3944  C Line 4078  C
4078        PARAMETER ( CH      = 2.5828674 )        PARAMETER ( CH      = 2.5828674 )
4079        PARAMETER ( CM      = CH / D1   )        PARAMETER ( CM      = CH / D1   )
4080    
4081          integer istnlv
4082        DIMENSION RI(irun,NLEV), STRT(irun,NLEV), DW2(irun,NLEV)        integer i
4083        DIMENSION XL(irun,NLEV),  ZKM(irun,NLEV), ZKH(irun,NLEV)  
       DIMENSION QE(irun,NLEV),  QQE(irun,NLEV)  
       INTEGER INTSTB(irun,nlev)  
       DIMENSION     EE(irun,nlay-1),  RF(irun,nlay-1)  
 C  
4084        ISTNLV = irun * NLEV        ISTNLV = irun * NLEV
4085  C  C
4086  C *********************************************************************  C *********************************************************************
# Line 4021  C    ZKH           -         HEAT TRANSP Line 4151  C    ZKH           -         HEAT TRANSP
4151  C    P3            -         PRODUCTION RATE OF TURBULENT KINETIC ENERG  C    P3            -         PRODUCTION RATE OF TURBULENT KINETIC ENERG
4152  C  C
4153  C**********************************************************************  C**********************************************************************
4154  C        implicit none
4155  C  
4156    C Argument list Declarations
4157          integer nlev,nlay,irun
4158          _RL Q(irun,NLEV),XL(irun,NLEV),STRT(irun,NLEV)
4159          _RL DW2(irun,NLEV)
4160          INTEGER INTSTB(irun,nlay), INTQ(irun,nlay)
4161          _RL ZKM(irun,NLEV),ZKH(irun,NLEV),P3(irun,NLEV)
4162    
4163    C Local Variables
4164          _RL a1,a2,a4,c1,a5,a3,b1,b2,b3,ff2,ff3,ff4
4165        PARAMETER ( A1      =   0.92    )        PARAMETER ( A1      =   0.92    )
4166        PARAMETER ( A2      =   0.74    )        PARAMETER ( A2      =   0.74    )
4167        PARAMETER ( A4      = 6. * A1 * A1)        PARAMETER ( A4      = 6. * A1 * A1)
# Line 4037  C Line 4175  C
4175        PARAMETER ( FF3     = (3.*A2*B2) - (9.*A2*A2 )  )        PARAMETER ( FF3     = (3.*A2*B2) - (9.*A2*A2 )  )
4176        PARAMETER ( FF4     = (3.*A2*B2) + (12.*A1*A2 )  )        PARAMETER ( FF4     = (3.*A2*B2) + (12.*A1*A2 )  )
4177    
4178          _RL F2(irun,nlay-1),F3(irun,nlay-1)
4179        DIMENSION   Q(irun,NLEV),  XL(irun,NLEV), STRT(irun,NLEV)        _RL F4(irun,nlay-1),XQ(irun,nlay-1)
4180        DIMENSION DW2(irun,NLEV)  
4181        DIMENSION ZKM(irun,NLEV), ZKH(irun,NLEV),   P3(irun,NLEV)        integer istnlv
4182  C        integer i
       DIMENSION F2(irun,nlay-1),  F3(irun,nlay-1)  
       DIMENSION F4(irun,nlay-1),  XQ(irun,nlay-1)  
       INTEGER INTSTB(irun,nlay), INTQ(irun,nlay)  
4183  C  C
4184        ISTNLV = irun * NLEV        ISTNLV = irun * NLEV
4185  C  C
# Line 4115  C    DXX1G         -         SOURCE TERM Line 4250  C    DXX1G         -         SOURCE TERM
4250  C    DXX1G         -         SOURCE TERM FOR XX2 AT GROUND  C    DXX1G         -         SOURCE TERM FOR XX2 AT GROUND
4251  C  C
4252  C**********************************************************************  C**********************************************************************
4253          implicit none
4254    
4255    C Argument List Declarations
4256          integer nlev,itype,irun
4257          _RL XX1(irun,NLEV+1),XX2(irun,NLEV+1)
4258          _RL RHOKDZ(irun,NLEV),FLXFAC(irun,NLEV+1)
4259          _RL DXX1G(irun),DXX2G(irun)
4260          _RL epsl
4261  C  C
4262  C        _RL AA(irun,nlev), BB(irun,nlev), CC(irun,nlev+1)
4263          integer istnlv,istnm1,nlevp1,istnlx
4264          integer i
       DIMENSION    XX1(irun,NLEV+1),    XX2(irun,NLEV+1)  
       DIMENSION RHOKDZ(irun,NLEV)  , FLXFAC(irun,NLEV+1)  
       DIMENSION  DXX1G(irun)    ,     DXX2G(irun)  
 C  
       DIMENSION AA(irun,nlev), BB(irun,nlev), CC(irun,nlev+1)  
4265  C  C
4266        ISTNLV = irun * NLEV        ISTNLV = irun * NLEV
4267        ISTNM1 = ISTNLV - irun        ISTNM1 = ISTNLV - irun
# Line 4182  C Line 4320  C
4320        RETURN        RETURN
4321        END        END
4322        SUBROUTINE VTRI0 ( A,B,C,F,Y,K,irun)        SUBROUTINE VTRI0 ( A,B,C,F,Y,K,irun)
4323        DIMENSION A(irun,K),B(irun,K),C(irun,K),Y(irun,K+1)        implicit none
4324        DIMENSION F(irun,K)  
4325          integer k,irun
4326          _RL A(irun,K),B(irun,K),C(irun,K),Y(irun,K+1)
4327          _RL F(irun,K)
4328    
4329          integer i,L,Lm1
4330  C  C
4331        DO 9000 I = 1,irun        DO 9000 I = 1,irun
4332         A(I,1) = 1. / A(I,1)         A(I,1) = 1. / A(I,1)
# Line 4208  C Line 4351  C
4351        END        END
4352  C  C
4353        SUBROUTINE VTRI1 ( A,B,Y,K,irun)        SUBROUTINE VTRI1 ( A,B,Y,K,irun)
4354        DIMENSION A(irun,K),B(irun,K),Y(irun,K+1)        implicit none
4355    
4356          integer k,irun
4357          _RL A(irun,K),B(irun,K),Y(irun,K+1)
4358    
4359          integer i,L
4360  C  C
4361        DO 200 L = K,1,-1        DO 200 L = K,1,-1
4362         DO 9000 I = 1,irun         DO 9000 I = 1,irun
# Line 4220  C Line 4368  C
4368        END        END
4369  C  C
4370        SUBROUTINE VTRI2 ( A,B,C,F,Y,K,irun)        SUBROUTINE VTRI2 ( A,B,C,F,Y,K,irun)
4371        DIMENSION A(irun,K),B(irun,K),C(irun,K),F(irun,K)        implicit none
4372        DIMENSION Y(irun,K+1)  
4373          integer k,irun
4374          _RL A(irun,K),B(irun,K),C(irun,K),F(irun,K)
4375          _RL Y(irun,K+1)
4376    
4377          integer i,L
4378  C  C
4379        DO 100 L = 2,K        DO 100 L = 2,K
4380         DO 9000 I = 1,irun         DO 9000 I = 1,irun
# Line 4282  C    DPSIH         -         D PSIH Line 4435  C    DPSIH         -         D PSIH
4435  C    BITRIB        -         BIT ARRAY - '1' WHERE RIB1 = 0  C    BITRIB        -         BIT ARRAY - '1' WHERE RIB1 = 0
4436  C  C
4437  C**********************************************************************  C**********************************************************************
4438  C        implicit none
4439  C  
4440    C Argument List Declarations
4441          integer nn,irun,itype
4442          _RL VRIB1(IRUN),VRIB2(IRUN)
4443          _RL VWS1(IRUN),VWS2(IRUN),VZ1(IRUN),VUSTAR(IRUN)
4444          integer IWATER(IRUN)
4445          _RL VAPSIM(IRUN),VAPSIHG(IRUN)
4446          _RL VPSIH(IRUN),VPSIG(IRUN),VX(IRUN)
4447          _RL VX0(IRUN),VY(IRUN),VY0(IRUN)
4448          LOGICAL LWATER
4449          _RL VDZETA(IRUN),VDZ0(IRUN),VDPSIM(IRUN)
4450          _RL VDPSIH(IRUN)
4451          integer INTRIB(IRUN)
4452          _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
4453          _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
4454          _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
4455          _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
4456          _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
4457          _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
4458          _RL VDPSIMC(irun),VDPSIHC(irun)
4459    
4460    C Local Variables
4461          _RL xx0max,prfac,xpfac,difsqt,ustz0s,h0byz0,usth0s
4462        PARAMETER ( XX0MAX  =   1.49821 )        PARAMETER ( XX0MAX  =   1.49821 )
4463        PARAMETER ( PRFAC  = 0.595864   )        PARAMETER ( PRFAC  = 0.595864   )
4464        PARAMETER ( XPFAC  = .55        )          PARAMETER ( XPFAC  = .55        )  
# Line 4292  C Line 4467  C
4467        PARAMETER ( H0BYZ0 =    30.0    )        PARAMETER ( H0BYZ0 =    30.0    )
4468        PARAMETER ( USTH0S =  H0BYZ0*USTZ0S )        PARAMETER ( USTH0S =  H0BYZ0*USTZ0S )
4469    
4470        DIMENSION VRIB1(IRUN),VRIB2(IRUN)        integer VINT1(irun),VINT2(irun)
4471        DIMENSION VWS1(IRUN),VWS2(IRUN),VZ1(IRUN),VUSTAR(IRUN)        _RL getcon,vk,bmdl,b2uhs
4472        DIMENSION VAPSIM(IRUN),VAPSIHG(IRUN)        integer i
       DIMENSION VPSIH(IRUN),VPSIG(IRUN),VX(IRUN)  
       DIMENSION VX0(IRUN),VY(IRUN),VY0(IRUN)  
       DIMENSION VDZETA(IRUN),VDZ0(IRUN),VDPSIM(IRUN)  
       DIMENSION VDPSIH(IRUN)  
       DIMENSION IWATER(IRUN),INTRIB(IRUN)  
       LOGICAL LWATER  
       DIMENSION VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)  
       DIMENSION VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)  
       DIMENSION VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)  
       DIMENSION VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)  
       DIMENSION VXNUM(irun),VDZETA1(irun),VDZETA2(irun)  
       DIMENSION VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)  
       DIMENSION VDPSIMC(irun),VDPSIHC(irun)  
 C  
   
   
       DIMENSION VINT1(irun),VINT2(irun)  
 C  
4473  C  C
4474        vk = getcon('VON KARMAN')        vk = getcon('VON KARMAN')
4475        BMDL    = VK * XPFAC * PRFAC / DIFSQT        BMDL    = VK * XPFAC * PRFAC / DIFSQT
# Line 4570  C        COMPUTE ROUGHNESS LENGTH FOR OC Line 4727  C        COMPUTE ROUGHNESS LENGTH FOR OC
4727  C          BASED ON FUNCTIONS OF LARGE AND POND  C          BASED ON FUNCTIONS OF LARGE AND POND
4728  C          AND OF KONDO --- DESIGNED FOR K = .4  C          AND OF KONDO --- DESIGNED FOR K = .4
4729  C *********************************************************************  C *********************************************************************
4730  C        implicit none
4731  C  
4732    C Argument List Delcarations
4733          integer irun
4734          _RL VZSEA(IRUN),VUSTAR(IRUN),VDZSEA(IRUN)
4735          integer IWATER(IRUN)
4736          LOGICAL LDZSEA
4737    
4738    C Local Variables
4739          _RL USTMX1,USTMX2,USTMX3
4740        PARAMETER ( USTMX1 =   1.14973  )        PARAMETER ( USTMX1 =   1.14973  )
4741        PARAMETER ( USTMX2 =   0.381844 )        PARAMETER ( USTMX2 =   0.381844 )
4742        PARAMETER ( USTMX3 =   0.0632456)        PARAMETER ( USTMX3 =   0.0632456)
4743    
4744        DIMENSION VZSEA(IRUN),VUSTAR(IRUN),VDZSEA(IRUN)        _RL AA(IRUN,5),TEMP(IRUN)
4745        DIMENSION IWATER(IRUN)        integer INT2(IRUN),INT3(IRUN),INT4(IRUN)
4746        DIMENSION AA(IRUN ,5),TEMP(IRUN)        integer i,k
4747        LOGICAL LDZSEA  
4748        DIMENSION INT2(IRUN ), INT3(IRUN ), INT4(IRUN)        _RL AA1(5),AA2(5),AA3(5),AA4(5)
 C  
       DIMENSION AA1(5),AA2(5),AA3(5),AA4(5)  
4749        DATA AA1/.2030325E-5,0.0,0.0,0.0,0.0/        DATA AA1/.2030325E-5,0.0,0.0,0.0,0.0/
4750        DATA AA2/-0.402451E-08,0.239597E-04,0.117484E-03,0.191918E-03,        DATA AA2/-0.402451E-08,0.239597E-04,0.117484E-03,0.191918E-03,
4751       1         0.395649E-04/       1         0.395649E-04/
# Line 4665  C Line 4828  C
4828  C  C
4829        RETURN        RETURN
4830        END        END
       subroutine pntquants_turb(nlaygcm)  
 C**********************************************************************  
 C  Subroutine to initialize the list of quantities for turbulence  
 C    package point by point diagnostic output  
 C**********************************************************************  
       implicit none  
       integer nlaygcm  
       character * 40 name  
       integer nlev  
   
       name = 'mid pressure'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'edge pressure'  
       nlev = nlaygcm+1  
       call pntquants(name,nlev)  
       name = 'mid p ** kappa'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'edge p ** kappa'  
       nlev = nlaygcm+1  
       call pntquants(name,nlev)  
       name = 'theta before'  
       nlev = nlaygcm+1  
       call pntquants(name,nlev)  
       name = 'theta virtual before'  
       nlev = nlaygcm+1  
       call pntquants(name,nlev)  
       name = 'q before'  
       nlev = nlaygcm+1  
       call pntquants(name,nlev)  
       name = 'u wind before'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'v wind before'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'heat cap ground'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'latent heat at ground'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'net surface rad'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'background heat transfer'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'tke before'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'richardson number before'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'theta after'  
       nlev = nlaygcm+1  
       call pntquants(name,nlev)  
       name = 'theta virtual after'  
       nlev = nlaygcm+1  
       call pntquants(name,nlev)  
       name = 'q after'  
       nlev = nlaygcm+1  
       call pntquants(name,nlev)  
       name = 'u wind after'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'v wind after'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'tke after'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'richardson number after'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'trb u flux'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'trb v flux'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'trb t flux'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'trb q flux'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'eddy coef mom'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'eddy coef heat'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'length scale'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'layer heights'  
       nlev = nlaygcm  
       call pntquants(name,nlev)  
       name = 'q ground'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'rib initial'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'cu initial'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'ct initial'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'ustar initial'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'rho sfc initial'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'z0 initial'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'zeta initial'  
       nlev = 1  
       call pntquants(name,nlev)  
       name = 'beta coeff'  
       nlev = 1  
       call pntquants(name,nlev)  
       return  
       end  
4831    
4832        subroutine seaice ( nocean, timstp, hice,        subroutine seaice ( nocean, timstp, hice,
4833       .                     eturb,  dedtc,         .                     eturb,  dedtc,  
# Line 4803  C*************************************** Line 4837  C***************************************
4837       .                       pke,  seaic, tc, qa )       .                       pke,  seaic, tc, qa )
4838        implicit none        implicit none
4839        integer  nocean        integer  nocean
4840        real timstp        _RL timstp
4841        real eturb(nocean),dedtc(nocean),hsturb(nocean),dhsdtc(nocean)        _RL eturb(nocean),dedtc(nocean),hsturb(nocean),dhsdtc(nocean)
4842        real swnet(nocean),lwnet(nocean),  dst4(nocean)        _RL swnet(nocean),lwnet(nocean),  dst4(nocean)
4843        real  qice(nocean),dqice(nocean)        _RL  qice(nocean),dqice(nocean)
4844        real   pke(nocean),   tc(nocean),    qa(nocean)        _RL   pke(nocean),   tc(nocean),    qa(nocean)
4845        real seaic(nocean)        _RL seaic(nocean)
4846    
4847  C  rho*C = 1.93e6 J/(m**3 K) ; Peixoto & Oort  C  rho*C = 1.93e6 J/(m**3 K) ; Peixoto & Oort
4848        real, parameter ::  rhoC = 1.93e6          _RL rhoC
4849          parameter (rhoC = 1.93e6)
4850    
4851        real faceps,getcon,latent,codt,deltg,hice        _RL faceps,getcon,latent,codt,deltg,hice
4852        integer i        integer i
4853    
4854        faceps = getcon('EPSFAC')        faceps = getcon('EPSFAC')
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22