/[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.13 by molod, Fri Jul 23 22:32:28 2004 UTC revision 1.47 by jmc, Wed Apr 1 23:58:22 2009 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "PACKAGES_CONFIG.h"  #include "FIZHI_OPTIONS.h"
 #include "CPP_OPTIONS.h"  
5        subroutine turbio (im,jm,nlay,istrip,nymd,nhms,bi,bj        subroutine turbio (im,jm,nlay,istrip,nymd,nhms,bi,bj
6       1 ,ndturb,ptop, pz, uz, vz, tz, qz, ntracers,ptracers       1 ,ndturb,nltop,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,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
# Line 14  C $Name$ Line 13  C $Name$
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,myid)       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
19  c  C
20  c  input:  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
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)
27  c      ndturb  - turbulence time step integer in HHMMSS format  C      ndturb  - turbulence time step integer in HHMMSS format
28  c      ptop    - model top pressure - rigid lid assumed - real  C      nltop   - Top level at which to allow turbulence
29  c      pz      - surface pressure minus ptop in mb - real[lon,lat]  C      ptop    - model top pressure - rigid lid assumed - real
30  c      uz      - zonal wind in m/sec - real[lon,lat,level]  C      pz      - surface pressure minus ptop in mb - real[lon,lat]
31  c      vz      - meridional wind in m/sec - real[lon,lat,level]  C      uz      - zonal wind in m/sec - real[lon,lat,level]
32  c      tz      - model theta (theta [deg K]/p0**k) - real[lon,lat,level]  C      vz      - meridional wind in m/sec - real[lon,lat,level]
33  c      qz      - specific humidity in kg/kg - real[lon,lat,level]  C      tz      - model theta (theta [deg K]/p0**k) - real[lon,lat,level]
34  c      ntracers- total number of tracers - integer  C      qz      - specific humidity in kg/kg - real[lon,lat,level]
35  c      ptracers- number of permanent tracers - integer  C      ntracers- total number of tracers - integer
36  c      pkht    - pressure[mb]**k at bottom edges of levels - real[lon,lat,level]  C      ptracers- number of permanent tracers - integer
37  c      fracland- not being used - real[lon,lat]  C      pkht    - pressure[mb]**k at bottom edges of levels - real[lon,lat,level]
38  c      landtype- not being used - integer[lon,lat]  C      fracland- not being used - real[lon,lat]
39  c      nchp    - nchplnd<nchp - total no chips (ocean too) - integer  C      landtype- not being used - integer[lon,lat]
40  c      nchplnd - <=nchp - number of land chips - integer  C      nchp    - nchplnd<nchp - total no chips (ocean too) - integer
41  c      chfr    - chip fraction - real[nchp]  C      nchplnd - <=nchp - number of land chips - integer
42  c      chlt    - tile space latitude array - real[nchp]  C      chfr    - chip fraction - real[nchp]
43  c      chlon   - tile space longitude array - real[nchp]  C      chlt    - tile space latitude array - real[nchp]
44  c      igrd    - tile space grid number - integer[nchp]  C      chlon   - tile space longitude array - real[nchp]
45  c      ityp    - tile space vegetation type - integer[nchp]  C      igrd    - tile space grid number - integer[nchp]
46  c      alai    - leaf area index - real[nchp]  C      ityp    - tile space vegetation type - integer[nchp]
47  c      agrn    - greenness fraction - real[nchp]  C      alai    - leaf area index - real[nchp]
48  c      thkz    - sea ice thickness in m (0. for no ice) - real[lon,lat]  C      agrn    - greenness fraction - real[nchp]
49  c      tprof   - logical flag for point by point diagnostic output  C      thkz    - sea ice thickness in m (0. for no ice) - real[lon,lat]
50  c      ndiagsiz- number of diagnostic 2-D arrays allocated  C      tprof   - logical flag for point by point diagnostic output
51  c      ndlsm   - number of tile diagnostic  C      ndiagsiz- number of diagnostic 2-D arrays allocated
52  c      radlwg  - net longwave flux at ground (up-down) in w/m**2 - real[lon,lat]  C      ndlsm   - number of tile diagnostic
53  c      st4     - upward longwave flux at ground in w/m**2 - real[lon,lat]  C      radlwg  - net longwave flux at ground (up-down) in w/m**2 - real[lon,lat]
54  c      dst4    - delta-sigma-T**4, ie, derivative of upward lw flux at  C      st4     - upward longwave flux at ground in w/m**2 - real[lon,lat]
55  c                ground with respect to ground Temperature - real[lon,lat]  C      dst4    - delta-sigma-T**4, ie, derivative of upward lw flux at
56  c      radswg  - net shortwave flux at ground (down-up) NON-DIM - real[lon,lat]  C                ground with respect to ground Temperature - real[lon,lat]
57  c                  {NOTE: this field is divided by the incident shortwave  C      radswg  - net shortwave flux at ground (down-up) NON-DIM - real[lon,lat]
58  c                     at the top of the atmosphere to non-dimensionalize]  C                  {NOTE: this field is divided by the incident shortwave
59  c      radswt  - incident shortwave at top of atmos in W/m**2 - real[lon,lat]  C                     at the top of the atmosphere to non-dimensionalize]
60  c      fdifpar - incident diffuse-beam PAR at surface in W/m**2 - real[lon,lat]  C      radswt  - incident shortwave at top of atmos in W/m**2 - real[lon,lat]
61  c      fdirpar - incident direct-beam PAR at surface in W/m**2 - real[lon,lat]  C      fdifpar - incident diffuse-beam PAR at surface in W/m**2 - real[lon,lat]
62  c      rainlsp - large-scale (frontal,supersat) rainfall in mm/sec - real[lon,lat]  C      fdirpar - incident direct-beam PAR at surface in W/m**2 - real[lon,lat]
63  c      rainconv- convective rainfall rate in mm/sec - real[lon,lat]  C      rainlsp - large-scale (frontal,supersat) rainfall in mm/sec - real[lon,lat]
64  c      snowfall- total snowfall rate in mm/sec - real[lon,lat]  C      rainconv- convective rainfall rate in mm/sec - real[lon,lat]
65  c updated:  C      snowfall- total snowfall rate in mm/sec - real[lon,lat]
66  c      tke     - turbulent k.e. in m**2/s**2 - real[tiles,levels]  C updated:
67  c      tgz     - surface skin temperature in deg K - real[lon,lat]  C      tke     - turbulent k.e. in m**2/s**2 - real[tiles,levels]
68  c      tcanopy - canopy temperature in deg K real[tiles]  C      tgz     - surface skin temperature in deg K - real[lon,lat]
69  c                 (sea surface temp over the ocean tiles)  C      tcanopy - canopy temperature in deg K real[tiles]
70  c      ecanopy - canopy vapor pressure in mb real[tiles]  C                 (sea surface temp over the ocean tiles)
71  c                 (qstar at tground over the sea ice and ocean tiles)  C      ecanopy - canopy vapor pressure in mb real[tiles]
72  c      tdeep   - deep soil temp in deg K real[tiles]  C                 (qstar at tground over the sea ice and ocean tiles)
73  c      swetshal- shallow level moisture field capacity fraction real[tiles]  C      tdeep   - deep soil temp in deg K real[tiles]
74  c      swetroot- root level moisture field capacity fraction real[tiles]  C      swetshal- shallow level moisture field capacity fraction real[tiles]
75  c      swetdeep- deep soil level moisture field capacity fraction real[tiles]  C      swetroot- root level moisture field capacity fraction real[tiles]
76  c      snodep  - depth of snow pack in cm liquid water equiv real[tiles]  C      swetdeep- deep soil level moisture field capacity fraction real[tiles]
77  c      capac   - leaf canopy water reservoir in cm real[tiles]  C      snodep  - depth of snow pack in cm liquid water equiv real[tiles]
78  c output:  C      capac   - leaf canopy water reservoir in cm real[tiles]
79  c      duturb  - change in zonal wind component due to turbulent processes  C output:
80  c                per unit time in m/sec**2 - real[lon,lat,levels]  C      duturb  - change in zonal wind component due to turbulent processes
81  c      dvturb  - change in meridional wind component due to turbulent processes  C                per unit time in m/sec**2 - real[lon,lat,levels]
82  c                per unit time in m/sec**2 - real[lon,lat,levels]  C      dvturb  - change in meridional wind component due to turbulent processes
83  c      dtturb  - change in (model theta*pi) due to turbulent processes  C                per unit time in m/sec**2 - real[lon,lat,levels]
84  c                per unit time  - real[lon,lat,levels] !! pi is pressure-ptop  C      dtturb  - change in (model theta*pi) due to turbulent processes
85  c      dqturb  - change in (specific humidity*pi) due to turbulent processes  C                per unit time  - real[lon,lat,levels] !! pi is pressure-ptop
86  c                per unit time  - real[lon,lat,levels] !! pi is pressure-ptop  C      dqturb  - change in (specific humidity*pi) due to turbulent processes
87  c      qliqavelw - Moist   Turbulence Liquid Water   for Longwave  - real[lon,lat,levels]  C                per unit time  - real[lon,lat,levels] !! pi is pressure-ptop
88  c      qliqavesw - Moist   Turbulence Liquid Water   for Shortwave - real[lon,lat,levels]  C      qliqavelw - Moist   Turbulence Liquid Water   for Longwave  - real[lon,lat,levels]
89  c      fccavelw  - Moist   Turbulence Cloud Fraction for Longwave  - real[lon,lat,levels]  C      qliqavesw - Moist   Turbulence Liquid Water   for Shortwave - real[lon,lat,levels]
90  c      fccavesw  - Moist   Turbulence Cloud Fraction for Shortwave - real[lon,lat,levels]  C      fccavelw  - Moist   Turbulence Cloud Fraction for Longwave  - real[lon,lat,levels]
91  c      qqgrid    - Gridded Turbulent Kinetic Energy                - real[lon,lat,levels]  C      fccavesw  - Moist   Turbulence Cloud Fraction for Shortwave - real[lon,lat,levels]
92  c-----------------------------------------------------------------------  C      qqgrid    - Gridded Turbulent Kinetic Energy                - real[lon,lat,levels]
93    C-----------------------------------------------------------------------
94        implicit none        implicit none
95    
96  #ifdef ALLOW_USE_MPI  #ifdef ALLOW_USE_MPI
97  #include "mpif.h"  #include "mpif.h"
98  #endif  #endif
99    
100  #ifdef ALLOW_DIAGNOSTICS        integer im,jm,nlay,istrip,nymd,nhms,bi,bj,ndturb,nltop
 #include "SIZE.h"  
 #include "diagnostics_SIZE.h"  
 #include "diagnostics.h"  
 #endif  
   
       integer im,jm,nlay,istrip,nymd,nhms,bi,bj,ndturb  
101        integer ntracers, ptracers        integer ntracers, ptracers
102        integer nchp,nchptot,nchplnd        integer nchp,nchptot,nchplnd
103        real ptop        _RL ptop
104        real pz(im,jm),uz(im,jm,nlay),vz(im,jm,nlay),tz(im,jm,nlay)  c     _RL pz(im,jm),uz(im,jm,nlay),vz(im,jm,nlay),tz(im,jm,nlay)
105        real qz(im,jm,nlay,ntracers)  c     _RL qz(im,jm,nlay,ntracers)
106        real plz(im,jm,nlay),plze(im,jm,nlay+1),dpres(im,jm,nlay)  c     _RL plz(im,jm,nlay),plze(im,jm,nlay+1),dpres(im,jm,nlay)
107        real pkht(im,jm,nlay+1),pkz(im,jm,nlay)  c     _RL pkht(im,jm,nlay+1),pkz(im,jm,nlay)
108        real ctmt(nchp),xxmt(nchp),yymt(nchp),zetamt(nchp)        _RL pz(im*jm,1),uz(im*jm,1,nlay),vz(im*jm,1,nlay)
109        real xlmt(nchp,nlay),khmt(nchp,nlay),tke(nchp,nlay)        _RL tz(im*jm,1,nlay),qz(im*jm,1,nlay,ntracers)
110        real tgz(im, jm),fracland(im,jm)        _RL plz(im*jm,1,nlay),plze(im*jm,1,nlay+1),dpres(im*jm,1,nlay)
111        integer landtype(im,jm)        _RL pkht(im*jm,1,nlay+1),pkz(im*jm,1,nlay)
112        real tcanopy(nchp),tdeep(nchp),ecanopy(nchp),swetshal(nchp)        _RL ctmt(nchp),xxmt(nchp),yymt(nchp),zetamt(nchp)
113        real swetroot(nchp),swetdeep(nchp),snodep(nchp),capac(nchp)        _RL xlmt(nchp,nlay),khmt(nchp,nlay),tke(nchp,nlay)
114        real chfr(nchp),chlt(nchp),chlon(nchp)  c     _RL tgz(im, jm),fracland(im,jm)
115    c     integer landtype(im,jm)
116          _RL tgz(im*jm,1),fracland(im*jm,1)
117          integer landtype(im*jm,1)
118          _RL tcanopy(nchp),tdeep(nchp),ecanopy(nchp),swetshal(nchp)
119          _RL swetroot(nchp),swetdeep(nchp),snodep(nchp),capac(nchp)
120          _RL chfr(nchp),chlt(nchp),chlon(nchp)
121        integer igrd(nchp),ityp(nchp)        integer igrd(nchp),ityp(nchp)
122        real alai(nchp),agrn(nchp),thkz(im,jm)  c     _RL alai(nchp),agrn(nchp),thkz(im,jm)
123          _RL alai(nchp),agrn(nchp),thkz(im*jm,1)
124        logical tprof        logical tprof
125        real duturb(im,jm,nlay),dvturb(im,jm,nlay)  c     _RL duturb(im,jm,nlay),dvturb(im,jm,nlay)
126        real dtturb(im,jm,nlay),dqturb(im,jm,nlay,ntracers)  c     _RL dtturb(im,jm,nlay),dqturb(im,jm,nlay,ntracers)
127        real st4(im,jm),dst4(im,jm)  c     _RL st4(im,jm),dst4(im,jm)
128        real radswg(im,jm),radswt(im,jm),radlwg(im,jm)  c     _RL radswg(im,jm),radswt(im,jm),radlwg(im,jm)
129        real fdifpar(im,jm),fdirpar(im,jm)  c     _RL fdifpar(im,jm),fdirpar(im,jm)
130        real rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)  c     _RL rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)
131        real tempref (im,jm)  c     _RL tempref (im,jm)
132          _RL duturb(im*jm,1,nlay),dvturb(im*jm,1,nlay)
133          _RL dtturb(im*jm,1,nlay),dqturb(im*jm,1,nlay,ntracers)
134          _RL st4(im*jm,1),dst4(im*jm,1)
135          _RL radswg(im*jm,1),radswt(im*jm,1),radlwg(im*jm,1)
136          _RL fdifpar(im*jm,1),fdirpar(im*jm,1)
137          _RL rainlsp(im*jm,1),rainconv(im*jm,1),snowfall(im*jm,1)
138          _RL tempref (im*jm,1)
139        integer imstturblw, imstturbsw        integer imstturblw, imstturbsw
140        real qliqavesw(im,jm,nlay),qliqavelw(im,jm,nlay)  c     _RL qliqavesw(im,jm,nlay),qliqavelw(im,jm,nlay)
141        real fccavelw (im,jm,nlay),fccavesw (im,jm,nlay)  c     _RL fccavelw (im,jm,nlay),fccavesw (im,jm,nlay)
142        real qqgrid   (im,jm,nlay)  c     _RL qqgrid   (im,jm,nlay)
143          _RL qliqavesw(im*jm,1,nlay),qliqavelw(im*jm,1,nlay)
144          _RL fccavelw (im*jm,1,nlay),fccavesw (im*jm,1,nlay)
145          _RL qqgrid   (im*jm,1,nlay)
146        integer myid        integer myid
147    
148  C Local Variables  C Local Variables
149    
150        integer numstrips        integer numstrips
151        integer ijall        integer ijall
152        real fmu,hice,tref,pref,cti,ed        _RL fmu,hice,tref,pref,cti,ed
153  C Set fmu and ed to zero for no background diffusion  C Set fmu and ed to zero for no background diffusion
154        parameter ( fmu    = 0.00000    )        parameter ( fmu    = 0.00000    )
155        parameter ( hice   =   300.     )        parameter ( hice   =   300.     )
156        parameter ( tref   =   258.     )          parameter ( tref   =   258.     )
157        parameter ( pref   =   500.     )        parameter ( pref   =   500.     )
158        parameter ( cti    =   0.0052   )        parameter ( cti    =   0.0052   )
159        parameter ( ed     =   0.0      )        parameter ( ed     =   0.0      )
160    
161        real qliqtmp(im,jm,nlay)  c     _RL qliqtmp(im,jm,nlay)
162        real  fcctmp(im,jm,nlay)  c     _RL  fcctmp(im,jm,nlay)
163        real tmpdiag(im,jm)  c     _RL tmpdiag(im,jm)
164        real   thtgz(im*jm)        _RL qliqtmp(im*jm,1,nlay)
165          _RL  fcctmp(im*jm,1,nlay)
166          _RL tmpdiag(im*jm,1), tmpFac
167          _RL   thtgz(im*jm)
168          logical  diagnostics_is_on
169          external diagnostics_is_on
170    
171        integer nland        integer nland
172        real alwcoeff(nchp),blwcoeff(nchp)        _RL alwcoeff(nchp),blwcoeff(nchp)
173        real netsw(nchp)        _RL netsw(nchp)
174        real cnvprec(nchp),lsprec(nchp)        _RL cnvprec(nchp),lsprec(nchp)
175        real snowprec(nchp)        _RL snowprec(nchp)
176        real pardiff(nchp),pardirct(nchp)        _RL pardiff(nchp),pardirct(nchp)
177        real pmsc(nchp)        _RL pmsc(nchp)
178        real netlw(nchp)        _RL netlw(nchp)
179        real sqscat(nchp), rsoil1(nchp)        _RL sqscat(nchp), rsoil1(nchp)
180        real rsoil2(nchp)        _RL rsoil2(nchp)
181        real rdc(nchp),u2fac(nchp)        _RL rdc(nchp),u2fac(nchp)
182        real z2ch(nchp)        _RL z2ch(nchp)
183        real zoch(nchp),cdrc(nchp)        _RL zoch(nchp),cdrc(nchp)
184        real cdsc(nchp)        _RL cdsc(nchp)
185        real dqsdt(nchp)        _RL dqsdt(nchp)
186        real tground(nchp),qground(nchp)        _RL tground(nchp),qground(nchp)
187        real utility(nchp)        _RL utility(nchp)
188        real    qice(nchp)        _RL    qice(nchp)
189        real   dqice(nchp)        _RL   dqice(nchp)
190    
191        real dumsc(nchp,nlay),dvmsc(nchp,nlay)        _RL dumsc(nchp,nlay),dvmsc(nchp,nlay)
192        real dtmsc(nchp,nlay),dqmsc(nchp,nlay,ntracers)        _RL dtmsc(nchp,nlay),dqmsc(nchp,nlay,ntracers)
193    
194        real shg(nchp),z0(nchp),icethk(nchp)        _RL shg(nchp),z0(nchp),icethk(nchp)
195        integer water(nchp)        integer water(nchp)
196    
197        real lats(istrip),lons(istrip),cosz(istrip),icest(istrip)        _RL lats(istrip),lons(istrip),cosz(istrip),icest(istrip)
198        real rainls(istrip),raincon(istrip),newsnow(istrip)        _RL rainls(istrip),raincon(istrip),newsnow(istrip)
199        real pardf(istrip),pardr(istrip),swnet(istrip)        _RL pardf(istrip),pardr(istrip),swnet(istrip)
200        real hlwdwn(istrip),alwrad(istrip),blwrad(istrip),tmpnlay(istrip)        _RL hlwdwn(istrip),alwrad(istrip),blwrad(istrip)
201        real laistrip(istrip),grnstrip(istrip),z2str(istrip),cd(istrip)        _RL tmpnlay(istrip)
202        real scatstr(istrip), rs1str(istrip), rs2str(istrip)        _RL laistrip(istrip),grnstrip(istrip),z2str(istrip),cd(istrip)
203        real rdcstr(istrip),u2fstr(istrip),dqsdtstr(istrip)        _RL scatstr(istrip), rs1str(istrip), rs2str(istrip)
204        real eturb(istrip),dedqa(istrip),dedtc(istrip)        _RL rdcstr(istrip),u2fstr(istrip),dqsdtstr(istrip)
205        real hsturb(istrip),dhsdqa(istrip),dhsdtc(istrip)        _RL eturb(istrip),dedqa(istrip),dedtc(istrip)
206        real savetc(istrip),saveqa(istrip),lwstrip(istrip)        _RL hsturb(istrip),dhsdqa(istrip),dhsdtc(istrip)
207        real chfrstr(istrip),psurf(istrip),shgstr(istrip)        _RL savetc(istrip),saveqa(istrip),lwstrip(istrip)
208          _RL chfrstr(istrip),psurf(istrip),shgstr(istrip)
209        integer types(istrip),igrdstr(istrip)        integer types(istrip),igrdstr(istrip)
210        real evap(istrip),shflux(istrip),runoff(istrip),bomb(istrip)        _RL evap(istrip),shflux(istrip),runoff(istrip),bomb(istrip)
211        real eint(istrip),esoi(istrip),eveg(istrip),esno(istrip)        _RL eint(istrip),esoi(istrip),eveg(istrip),esno(istrip)
212        real smelt(istrip),hlatn(istrip),hlwup(istrip),gdrain(istrip)        _RL smelt(istrip),hlatn(istrip),hlwup(istrip),gdrain(istrip)
213        real runsrf(istrip),fwsoil(istrip),evpot(istrip)        _RL runsrf(istrip),fwsoil(istrip),evpot(istrip)
214        real strdg1(istrip),strdg2(istrip),strdg3(istrip),strdg4(istrip)        _RL strdg1(istrip),strdg2(istrip),strdg3(istrip),strdg4(istrip)
215        real strdg5(istrip),strdg6(istrip),strdg7(istrip),strdg8(istrip)        _RL strdg5(istrip),strdg6(istrip),strdg7(istrip),strdg8(istrip)
216        real strdg9(istrip),tmpstrip(istrip),qicestr(istrip)        _RL strdg9(istrip),tmpstrip(istrip),qicestr(istrip)
217        real dqicestr(istrip)        _RL dqicestr(istrip)
218    
219        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)
220        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)
221        real tracers(istrip,nlay+1,ntracers)        _RL tracers(istrip,nlay+1,ntracers)
222        real dpstr(istrip,nlay),pke(istrip,nlay+1)        _RL dpstr(istrip,nlay),pke(istrip,nlay+1)
223        real pk(istrip,nlay), qq(istrip,nlay),   p(istrip,nlay)        _RL pk(istrip,nlay), qq(istrip,nlay),   p(istrip,nlay)
224        real sri(istrip,nlay), skh(istrip,nlay), skm(istrip,nlay)        _RL sri(istrip,nlay), skh(istrip,nlay), skm(istrip,nlay)
225        real stuflux(istrip,nlay), stvflux(istrip,nlay)        _RL stuflux(istrip,nlay), stvflux(istrip,nlay)
226        real sttflux(istrip,nlay), stqflux(istrip,nlay)        _RL sttflux(istrip,nlay), stqflux(istrip,nlay)
227        real frqtrb(istrip,nlay-1)        _RL frqtrb(istrip,nlay-1)
228        real dshdthg(istrip,nlay),dthdthg(istrip,nlay)        _RL dshdthg(istrip,nlay),dthdthg(istrip,nlay)
229        real dshdshg(istrip,nlay),dthdshg(istrip,nlay)        _RL dshdshg(istrip,nlay),dthdshg(istrip,nlay)
230        real transth(istrip,nlay), transsh(istrip,nlay)        _RL transth(istrip,nlay), transsh(istrip,nlay)
231    
232        real tc(istrip),td(istrip),qa(istrip)        _RL tc(istrip),td(istrip),qa(istrip)
233        real swet1(istrip),swet2(istrip),swet3(istrip)        _RL swet1(istrip),swet2(istrip),swet3(istrip)
234        real capacity(istrip),snowdepth(istrip)        _RL capacity(istrip),snowdepth(istrip)
235        real stz0(istrip)        _RL stz0(istrip)
236        real stdiag(istrip)        _RL stdiag(istrip)
237        real tends(istrip),sustar(istrip), sz0(istrip),pbldpth(istrip)        _RL tends(istrip),sustar(istrip), sz0(istrip),pbldpth(istrip)
238        real sct(istrip), scu(istrip), swinds(istrip)        _RL sct(istrip), scu(istrip), swinds(istrip)
239        real stu2m(istrip),stv2m(istrip),stt2m(istrip),stq2m(istrip)        _RL stu2m(istrip),stv2m(istrip),stt2m(istrip),stq2m(istrip)
240        real stu10m(istrip),stv10m(istrip),stt10m(istrip),stq10m(istrip)        _RL stu10m(istrip),stv10m(istrip),stt10m(istrip),stq10m(istrip)
241        integer  stwatr(istrip)        integer  stwatr(istrip)
242        real  wspeed(istrip)        _RL  wspeed(istrip)
243    
244        real ctsave(istrip),xxsave(istrip),yysave(istrip),zetasave(istrip)        _RL ctsave(istrip),xxsave(istrip),yysave(istrip)
245        real xlsave(istrip,nlay),khsave(istrip,nlay)        _RL zetasave(istrip)
246        real qliq(istrip,nlay),turbfcc(istrip,nlay)        _RL xlsave(istrip,nlay),khsave(istrip,nlay)
247        real qliqmsc(nchp,nlay),fccmsc(nchp,nlay)        _RL qliq(istrip,nlay),turbfcc(istrip,nlay)
248          _RL qliqmsc(nchp,nlay),fccmsc(nchp,nlay)
249    
250        integer ndlsm        integer ndlsm
251        parameter ( ndlsm = 40)  c     parameter ( ndlsm = 1)
252        real qdiaglsm(nchp,ndlsm)        parameter ( ndlsm = 38)
253          _RL qdiaglsm(nchp,ndlsm)
254    
255          _RL pi,secday,sdayopi2,rgas,akap,cp,alhl
256          _RL faceps,grav,caltoj,virtcon,getcon
257          _RL heatw,undef,timstp,delttrb,dttrb,ra
258          _RL edle,rmu,cltj10,atimstp,tice,const
259          integer istnp1,istnlay,itrtrb,i,L,nn,nt
260    c     integer j
261          integer nocean, nice
262          integer ndmoist,time_left, ndum0,ndum1
263          integer ntracedim
264          _RL    dtfac,timstp2,sum0
265    C logical begin flag - set to true to indicate a cold start
266          logical qbeg
267    
268        integer n,nsecf,nmonf,ndayf        integer n,nsecf,nmonf,ndayf
269        nsecf(n) = n/10000*3600 + mod(n,10000)/100* 60 + mod(n,100)        nsecf(n) = n/10000*3600 + mod(n,10000)/100* 60 + mod(n,100)
270        nmonf(n) = mod(n,10000)/100        nmonf(n) = mod(n,10000)/100
271        ndayf(n) = mod(n,100)        ndayf(n) = mod(n,100)
272    
       real pi,secday,sdayopi2,rgas,akap,cp,alhl  
       real faceps,grav,caltoj,virtcon,getcon  
       real heatw,undef,timstp,delttrb,dttrb,ra  
       real edle,rmu,cltj10,atimstp,tice,const  
       integer istnp1,istnlay,itrtrb,i,j,L,nn,nt  
       integer nocean, nice  
       integer ndmoist,time_left,ndum  
       integer ntracedim  
       real    dtfac,timstp2,sum  
 C logical begin flag - set to true to indicate a cold start  
       logical qbeg      
   
273  #ifdef CRAY  #ifdef CRAY
274  #ifdef f77  #ifdef f77
275  cfpp$ expand (qsat)  cfpp$ expand (qsat)
276  #endif  #endif
277  #endif  #endif
278    
279  c   compute variables that do not change  C   compute variables that do not change
280  c  
281         pi = 4.*atan(1.)         pi = 4.*atan(1.)
282         secday   = getcon('SDAY')         secday   = getcon('SDAY')
283         sdayopi2 = getcon('SDAY') / (pi*2.)         sdayopi2 = getcon('SDAY') / (pi*2.)
# Line 273  c Line 293  c
293         undef    = getcon('UNDEF')         undef    = getcon('UNDEF')
294         ntracedim= max(ntracers-ptracers,1)         ntracedim= max(ntracers-ptracers,1)
295    
296         call get_alarm ( 'moist',ndum,ndum,ndmoist,time_left )         call get_alarm ( 'moist',ndum0,ndum1,ndmoist,time_left )
297         timstp   = nsecf(ndturb)         timstp   = nsecf(ndturb)
298         timstp2  = nsecf(ndmoist)         timstp2  = nsecf(ndmoist)
299         dtfac    = min( 1.0, timstp/timstp2 )         dtfac    = min( 1.0 _d 0, timstp/timstp2 )
300    
301  c delttrb is the internal turbulence time step  C delttrb is the internal turbulence time step
302  c a value equal to ndturb means one internal iteration  C a value equal to ndturb means one internal iteration
303         delttrb = nsecf(ndturb)         delttrb = nsecf(ndturb)
304    
305         ijall    =   im * jm         ijall    =   im * jm
# Line 288  c a value equal to ndturb means one inte Line 308  c a value equal to ndturb means one inte
308         itrtrb   = ( timstp / delttrb ) + 0.1         itrtrb   = ( timstp / delttrb ) + 0.1
309         dttrb    =   timstp / float(itrtrb)         dttrb    =   timstp / float(itrtrb)
310         edle     =   ed * 0.2         edle     =   ed * 0.2
311    
312  c   coefficient of viscosity (background momentum diffusion)  C   coefficient of viscosity (background momentum diffusion)
313  c  C
314         rmu     = fmu * tref * rgas / pref         rmu     = fmu * tref * rgas / pref
315         cltj10  = 10. * caltoj         cltj10  = 10. * caltoj
316         atimstp = 1. / timstp         atimstp = 1. / timstp
317         tice = getcon('FREEZING-POINT')         tice = getcon('FREEZING-POINT')
318    
319  c **********************************************************************  C **********************************************************************
320  c   Check for Cold Start (if QQ is zero everywhere)  C   Check for Cold Start (if QQ is zero everywhere)
321  c **********************************************************************  C **********************************************************************
322          
323        qbeg = .false.        qbeg = .false.
324    
325        sum = 0.0        sum0 = 0.0
326        do L=1,nlay        do L=1,nlay
327        do n=1,nchptot        do n=1,nchptot
328        sum = sum + tke(n,L)        sum0 = sum0 + tke(n,L)
329        enddo        enddo
330        enddo        enddo
331    
332  #ifdef ALLOW_USE_MPI  #ifdef ALLOW_USE_MPI
333        call mpi_allreduce(sum,const,1,mpi_double_precision,mpi_sum,        call mpi_allreduce(sum0,const,1,mpi_double_precision,mpi_sum,
334       .                                                mpi_comm_world,n)       .                                                mpi_comm_world,n)
335  #else  #else
336        const = sum        const = sum0
337  #endif  #endif
338    
339        if( const.eq.0.0 ) then        if( const.eq.0.0 ) then
340        qbeg = .true.        qbeg = .true.
341            if( myid.eq.0 ) then            if( myid.eq.1 .and. bi.eq.1 ) then
342            print *            print *
343            print *, 'Warning!'            print *, 'Warning!'
344            print *, 'Turbulent Kinetic Energy has not been initialized.'            print *, 'Turbulent Kinetic Energy has not been initialized.'
# Line 327  c ************************************** Line 347  c **************************************
347            endif            endif
348        endif        endif
349    
350  c **********************************************************************  C **********************************************************************
351  c                            Initialization  C                            Initialization
352  c **********************************************************************  C **********************************************************************
         
 c Initialize diagnostic for ground temperature change  
 c ---------------------------------------------------  
       if(idtg.gt.0) then  
       do i =1,ijall  
       qdiag(i,1,idtg,bi,bj) = qdiag(i,1,idtg,bi,bj) - tgz(i,1)  
       enddo  
       endif  
353    
354  c **********************************************************************  C Initialize diagnostic for ground temperature change
355  c  entire turbulence and land surface package will run in 'tile space'    C ---------------------------------------------------
356  c       do conversion of model state variables to tile space  c     if(diagnostics_is_on('DTG     ',myid) ) then
357  c        (ocean points appended to tile space land point arrays)  c      do j =1,jm
358  c **********************************************************************  c      do i =1,im
359    c       tmpdiag(i,j) =  -tgz(i,j)
360    c      enddo
361    c      enddo
362    c      call diagnostics_fill(tmpdiag,'DTG     ',0,1,-3,bi,bj,myid)
363    c     endif
364          tmpFac = -1.
365          CALL DIAGNOSTICS_SCALE_FILL( tgz,tmpFac,1,
366         &                             'DTG     ',0,1,-3,bi,bj,myid)
367    
368    C **********************************************************************
369    C  entire turbulence and land surface package will run in 'tile space'
370    C       do conversion of model state variables to tile space
371    C        (ocean points appended to tile space land point arrays)
372    C **********************************************************************
373    
374        numstrips   = ((nchptot-1)/istrip) + 1        numstrips   = ((nchptot-1)/istrip) + 1
375    
376        call grd2msc(pz(1,1),im,jm,igrd,pmsc,nchp,nchptot)        call grd2msc(pz(1,1),im,jm,igrd,pmsc,nchp,nchptot)
377    
378        call grd2msc(tgz,im,jm,igrd,tground,nchp,nchptot)        call grd2msc(tgz,im,jm,igrd,tground,nchp,nchptot)
379        do i = 1,ijall        do i = 1,ijall
380         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))
381       1                         - dst4(i,1)* tgz(i,1)       1                         - dst4(i,1)* tgz(i,1)
382        enddo        enddo
383        call grd2msc(tmpdiag,im,jm,igrd,alwcoeff,nchp,nchptot)        call grd2msc(tmpdiag,im,jm,igrd,alwcoeff,nchp,nchptot)
# Line 384  C Call chpprm to get non-varying vegetat Line 411  C Call chpprm to get non-varying vegetat
411        call chpprm(nymd,nhms,nchp,nchplnd,chlt,ityp,alai,        call chpprm(nymd,nhms,nchp,nchplnd,chlt,ityp,alai,
412       1       agrn,zoch,z2ch,cdrc,cdsc,sqscat,u2fac,rsoil1,rsoil2,rdc)       1       agrn,zoch,z2ch,cdrc,cdsc,sqscat,u2fac,rsoil1,rsoil2,rdc)
413    
414  c **********************************************************************  C **********************************************************************
415  c ****                surface specification                         ****  C ****                surface specification                         ****
416  c **********************************************************************  C **********************************************************************
417    
418  c   set water  C   set water
419    
420        do i = 1,nchptot        do i = 1,nchptot
421         water(i) = 0         water(i) = 0
422         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
423        enddo        enddo
424    
425  c   roughness length z0  C   roughness length z0
426  c  C
427        do i =1,nchptot        do i =1,nchptot
428         if (icethk(i).gt.0.) then         if (icethk(i).gt.0.) then
429          z0(i) = 1.e-4          z0(i) = 1.e-4
# Line 407  c Line 434  c
434         endif         endif
435        enddo        enddo
436    
437  c Fill Array Tground with canopy temperatures over land tiles  C Fill Array Tground with canopy temperatures over land tiles
438  c  (it has sst from the tgz array over the sea ice and ocean tiles)  C  (it has sst from the tgz array over the sea ice and ocean tiles)
439    
440        do i = 1,nchplnd        do i = 1,nchplnd
441         tground(i) = tcanopy(i)         tground(i) = tcanopy(i)
# Line 421  C --------------------- Line 448  C ---------------------
448        call qsat ( tground(i),utility(i),shg(i),dqsdt(i),.true. )        call qsat ( tground(i),utility(i),shg(i),dqsdt(i),.true. )
449        enddo        enddo
450    
451  c Fill Array Qground with canopy air specific humidity over land tiles  C Fill Array Qground with canopy air specific humidity over land tiles
452  c  (it has qstar at tground over the sea ice and ocean tiles)  C  (it has qstar at tground over the sea ice and ocean tiles)
453    
454        do i = 1,nchplnd        do i = 1,nchplnd
455         qground(i) = ecanopy(i)         qground(i) = ecanopy(i)
# Line 431  c  (it has qstar at tground over the sea Line 458  c  (it has qstar at tground over the sea
458         qground(i) = shg(i)         qground(i) = shg(i)
459        enddo        enddo
460    
461  c Fill Array Swetshal with Value 1. over oceans and sea ice  C Fill Array Swetshal with Value 1. over oceans and sea ice
462        do i = nchplnd+1,nchptot        do i = nchplnd+1,nchptot
463         swetshal(i) = 1.         swetshal(i) = 1.
464        enddo        enddo
465    
466  c compute heat conduction through ice  C compute heat conduction through ice
467  c -----------------------------------  C -----------------------------------
468        const = ( cti / hice ) * cltj10        const = ( cti / hice ) * cltj10
469        do i =1,nchptot        do i =1,nchptot
470               qice(i) =  0.0               qice(i) =  0.0
# Line 448  c ----------------------------------- Line 475  c -----------------------------------
475         endif         endif
476        enddo        enddo
477    
478        if( iqice.gt.0 ) then        if(diagnostics_is_on('QICE    ',myid) ) then
479        do i =1,ijall         do i =1,ijall
480        tmpdiag(i,1) = 0.0          tmpdiag(i,1) = 0.0
481        enddo         enddo
482        call msc2grd (igrd,chfr,qice,nchp,nchptot,fracland,tmpdiag,im,jm)         call msc2grd (igrd,chfr,qice,nchp,nchptot,fracland,tmpdiag,im,jm)
483        do i =1,ijall         call diagnostics_fill(tmpdiag,'QICE    ',0,1,3,bi,bj,myid)
       qdiag(i,1,iqice,bi,bj) = qdiag(i,1,iqice,bi,bj) + tmpdiag(i,1)  
       enddo  
       nqice = nqice + 1  
484        endif        endif
485    
486  c**********************************************************************  C***********************************************************************
487  c                        loop over regions  C                        loop over regions
488  c**********************************************************************  C***********************************************************************
489    
490        do 2000 nn = 1, numstrips        do 2000 nn = 1, numstrips
491    
# Line 474  c*************************************** Line 498  c***************************************
498         call strip2tile(plze,igrd,pe,nchp,ijall,istrip,nlay+1,nn)         call strip2tile(plze,igrd,pe,nchp,ijall,istrip,nlay+1,nn)
499         call strip2tile(pkz,igrd,pk,nchp,ijall,istrip,nlay,nn)         call strip2tile(pkz,igrd,pk,nchp,ijall,istrip,nlay,nn)
500         call strip2tile(pkht,igrd,pke,nchp,ijall,istrip,nlay+1,nn)         call strip2tile(pkht,igrd,pke,nchp,ijall,istrip,nlay+1,nn)
501         do nt = 1,ntracers-ptracers  c      do nt = 1,ntracers-ptracers
502         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,
503       1                                             ijall,istrip,nlay,nn)  c    1                                             ijall,istrip,nlay,nn)
504         enddo  c      enddo
505    
506         call stripit  (z0,stz0,nchptot,nchp,istrip,1,nn)         call stripit  (z0,stz0,nchptot,nchp,istrip,1,nn)
507         call stripit  (tground,th(1,nlay+1),nchptot,nchp,istrip,1,nn)         call stripit  (tground,th(1,nlay+1),nchptot,nchp,istrip,1,nn)
# Line 528  c*************************************** Line 552  c***************************************
552         call stripit  (snodep,snowdepth,nchptot,nchp,istrip,1,nn)         call stripit  (snodep,snowdepth,nchptot,nchp,istrip,1,nn)
553         call stripit  (capac,capacity,nchptot,nchp,istrip,1,nn)         call stripit  (capac,capacity,nchptot,nchp,istrip,1,nn)
554    
555    #ifdef FIZHI_USE_FIXED_DAY
556           call astro ( 20040321,nhms,lats,lons,istrip,cosz,ra )
557    #else
558         call astro ( nymd,nhms,lats,lons,istrip,cosz,ra )         call astro ( nymd,nhms,lats,lons,istrip,cosz,ra )
559    #endif
560    
561  c we need to count up the land, sea ice and ocean points  C we need to count up the land, sea ice and ocean points
562        nocean = 0        nocean = 0
563        nland  = 0        nland  = 0
564        nice   = 0        nice   = 0
# Line 540  c we need to count up the land, sea ice Line 568  c we need to count up the land, sea ice
568         if( types(i).eq.100 .and. icest(i).gt.0.0 ) nice = nice + 1         if( types(i).eq.100 .and. icest(i).gt.0.0 ) nice = nice + 1
569        enddo        enddo
570    
571  c Disable following ISTRIP check for MPI version  C Zero out velocities at the bottom edge of the model
572  c ----------------------------------------------  C ---------------------------------------------------
573  c     if( (nland+nocean).ne.istrip ) then        do i =1,istrip
574  c     print *         u(i,nlay+1) = 0.
575  c     print *,'Error!'         v(i,nlay+1) = 0.
576  c     print *,'Problem Stripping Land/Ocean/Ice points in Turbulence'        enddo
 c     print *  
 c     stop  
 c     endif  
577    
578  c convert temperature of level nlay+1 to theta & value of sh at ground  C convert temperature of level nlay+1 to theta & value of sh at ground
579  c --------------------------------------------------------------------  C --------------------------------------------------------------------
580        do i =1,istrip        do i =1,istrip
581        th(i,nlay+1) = th(i,nlay+1) / pke(i,nlay+1)        th(i,nlay+1) = th(i,nlay+1) / pke(i,nlay+1)
582        sh(i,nlay+1) = qa(i)        sh(i,nlay+1) = qa(i)
583        enddo        enddo
584    
585        if(iqg.gt.0) then        if(diagnostics_is_on('QG      ',myid) ) then
586        do i=1,istrip        do i=1,istrip
587        tmpstrip(i) = sh(i,nlay+1)*1000        tmpstrip(i) = sh(i,nlay+1)*1000
588        enddo        enddo
589        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
590       1                        qdiag(1,1,iqg,bi,bj),ijall,1,nn,.false.)       . .false., 'QG      ', 1, 1, bi, bj, myid)
591        endif        endif
592    
593  c value of tracers at the ground  C value of tracers at the ground
594  c ------------------------------  C ------------------------------
595        do nt = 1,ntracers-ptracers  c     do nt = 1,ntracers-ptracers
596         do i = 1,istrip  c      do i = 1,istrip
597          tracers(i,nlay+1,nt) = 0.  c       tracers(i,nlay+1,nt) = 0.
598         enddo  c      enddo
599        enddo  c     enddo
600    
601  c compute virtual potential temperatures  C compute virtual potential temperatures
602  c --------------------------------------  C --------------------------------------
603        do L = 1,nlay+1        do L = 1,nlay+1
604        do i =1,istrip        do i =1,istrip
605        thv(i,L) = 1. + virtcon * sh(i,L)        thv(i,L) = 1. + virtcon * sh(i,L)
# Line 585  c -------------------------------------- Line 610  c --------------------------------------
610        sh(i,nlay+1) = qa(i)        sh(i,nlay+1) = qa(i)
611        enddo        enddo
612    
613  c zero out arrays for output of qliq and fcc  C zero out arrays for output of qliq and fcc
614        do L =1,nlay        do L =1,nlay
615        do i =1,istrip        do i =1,istrip
616        qliq(i,L) = 0.        qliq(i,L) = 0.
# Line 593  c zero out arrays for output of qliq and Line 618  c zero out arrays for output of qliq and
618        enddo        enddo
619        enddo        enddo
620    
621  c zero out fluxes and derivatives  C zero out fluxes and derivatives
622  c -------------------------------  C -------------------------------
623        do i = 1,istrip        do i = 1,istrip
624         eturb(i) = 0.         eturb(i) = 0.
625         scu(i) = 0.         scu(i) = 0.
# Line 605  c ------------------------------- Line 630  c -------------------------------
630         dhsdtc(i) = 0.         dhsdtc(i) = 0.
631        enddo        enddo
632    
633  c increment diagnostic arrays for quantities calculated before trbfl  C increment diagnostic arrays for quantities calculated before trbfl
634  c ------------------------------------------------------------------  C ------------------------------------------------------------------
635        do i =1,istrip        if(diagnostics_is_on('DTSRF   ',myid) ) then
636        stdiag(i) = ( thv(i,nlay+1)-thv(i,nlay) ) / pke(i,nlay+1)         do i=1,istrip
637        enddo          stdiag(i) = ( thv(i,nlay+1)-thv(i,nlay) ) / pke(i,nlay+1)
638        if(idtsrf.gt.0) then         enddo
639        call paste2grd(stdiag,igrd,chfrstr,istrip,nchp,         call diag_vegtile_fill (stdiag,igrd,chfrstr,istrip,nchp,nn,
640       1                      qdiag(1,1,idtsrf,bi,bj),ijall,1,nn,.false.)       . .false., 'DTSRF   ', 1, 1, bi, bj, myid)
641        endif        endif
642    
643  c call trbflx  C call trbflx
644  c -----------  C -----------
645        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,
646       1 tracers,ntracers-ptracers,ntracedim,dttrb,itrtrb,rmu,edle,qbeg,       1 tracers,ntracers-ptracers,ntracedim,dttrb,itrtrb,rmu,edle,qbeg,
647       2 tprof,stuflux,stvflux,sri,skh,skm,swinds,sustar,sz0,frqtrb,       2 tprof,stuflux,stvflux,sri,skh,skm,swinds,sustar,sz0,frqtrb,
648       3 pbldpth,sct,scu,stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,       3 pbldpth,sct,scu,stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,
649       4 stq10m,istrip,nlay,nymd,nhms,grav,cp,rgas,faceps,virtcon,undef,       4 stq10m,istrip,nlay,nltop,nymd,nhms,grav,cp,rgas,faceps,virtcon,
650       5 dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,       5 undef,dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,
651       6 hsturb,dhsdqa,dhsdtc,transth,transsh,       6 hsturb,dhsdqa,dhsdtc,transth,transsh,
652       7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)       7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)
653    
# Line 637  c ----------- Line 662  c -----------
662        call pastit (qliq  ,qliqmsc,istrip,nchp,nchptot,nlay,nn)        call pastit (qliq  ,qliqmsc,istrip,nchp,nchptot,nlay,nn)
663        call pastit (turbfcc,fccmsc,istrip,nchp,nchptot,nlay,nn)        call pastit (turbfcc,fccmsc,istrip,nchp,nchptot,nlay,nn)
664    
665  c  New diagnostic: potential evapotranspiration  C  New diagnostic: potential evapotranspiration
666        do i = 1,istrip        do i = 1,istrip
667         evpot(i) = transsh(i,nlay) * (shgstr(i) - sh(i,nlay))         evpot(i) = transsh(i,nlay) * (shgstr(i) - sh(i,nlay))
668        enddo        enddo
# Line 645  c  New diagnostic: potential evapotransp Line 670  c  New diagnostic: potential evapotransp
670  C**********************************************************************  C**********************************************************************
671  C   Call Land Surface Module  C   Call Land Surface Module
672  C**********************************************************************  C**********************************************************************
673    
674        do i = 1,istrip        do i = 1,istrip
675         savetc(i) = tc(i)         savetc(i) = tc(i)
676         saveqa(i) = qa(i)         saveqa(i) = qa(i)
677        enddo        enddo
678        do i = 1,istrip        do i = 1,istrip
679         cosz(i) = max(cosz(i),0.0001)         cosz(i) = max(cosz(i),0.0001 _d 0)
680         cd(i) = scu(i)*scu(i)         cd(i) = scu(i)*scu(i)
681         tmpnlay(i) = th(i,nlay)*pk(i,nlay)         tmpnlay(i) = th(i,nlay)*pk(i,nlay)
682         hlwdwn(i) = alwrad(i)+blwrad(i)*tc(i)-lwstrip(i)         hlwdwn(i) = alwrad(i)+blwrad(i)*tc(i)-lwstrip(i)
683         psurf(i) = pe(i,nlay+1)         psurf(i) = pe(i,nlay+1)
684         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))
685           if(wspeed(i) .lt. 1.e-20) wspeed(i) = 1.e-20
686  C Note: This LSM precip bug needs to be cleaned up  C Note: This LSM precip bug needs to be cleaned up
687  ccc   newsnow(i) = newsnow(i)*dtfac    ccc   newsnow(i) = newsnow(i)*dtfac
688  ccc   raincon(i) = raincon(i)*dtfac    ccc   raincon(i) = raincon(i)*dtfac
689  ccc   rainls (i) = rainls (i)*dtfac    ccc   rainls (i) = rainls (i)*dtfac
690        enddo        enddo
691    
692        do i = 1,istrip        do i = 1,istrip
# Line 693  ccc   rainls (i) = rainls (i)*dtfac Line 719  ccc   rainls (i) = rainls (i)*dtfac
719         fwsoil(i) = 0.         fwsoil(i) = 0.
720        enddo        enddo
721    
722  c**********************************************************************  C**********************************************************************
723  c   diagnostics: fill arrays for lsm input fields  C   diagnostics: fill arrays for lsm input fields
724  c**********************************************************************  C**********************************************************************
725        if(isnowfall.gt.0) then        if(diagnostics_is_on('SNOWFALL',myid) ) then
726        do i = 1,istrip        do i = 1,istrip
727        tmpstrip(i) = newsnow(i)*86400          tmpstrip(i) = newsnow(i)*86400
728        enddo        enddo
729        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
730       1                   qdiag(1,1,isnowfall,bi,bj),ijall,1,nn,.false.)       . .false., 'SNOWFALL', 1, 1, bi, bj, myid)
731        endif        endif
732        if(iraincon.gt.0) then        if(diagnostics_is_on('RAINCON ',myid) ) then
733        do i = 1,istrip        do i = 1,istrip
734        tmpstrip(i) = raincon(i)*86400          tmpstrip(i) = raincon(i)*86400
735        enddo        enddo
736        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
737       1                   qdiag(1,1,iraincon,bi,bj),ijall,1,nn,.false.)       . .false., 'RAINCON ', 1, 1, bi, bj, myid)
738        endif        endif
739        if(irainlsp.gt.0) then        if(diagnostics_is_on('RAINLSP ',myid) ) then
740        do i = 1,istrip        do i = 1,istrip
741        tmpstrip(i) = rainls(i)*86400          tmpstrip(i) = rainls(i)*86400
742        enddo        enddo
743        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
744       1                  qdiag(1,1,irainlsp,bi,bj),ijall,1,nn,.false.)       . .false., 'RAINLSP ', 1, 1, bi, bj, myid)
       endif  
       if(igreen.gt.0) then  
       call paste2grd(grnstrip,igrd,chfrstr,istrip,nchp,  
      1                    qdiag(1,1,igreen,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(ilai.gt.0) then  
       call paste2grd(laistrip,igrd,chfrstr,istrip,nchp,  
      1                    qdiag(1,1,ilai,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(ipardr.gt.0) then  
       call paste2grd(pardr,igrd,chfrstr,istrip,nchp,  
      1                 qdiag(1,1,ipardr,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(ipardf.gt.0) then  
       call paste2grd(pardf,igrd,chfrstr,istrip,nchp,  
      1                 qdiag(1,1,ipardf,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(idlwdtc.gt.0) then  
       call paste2grd(blwrad,igrd,chfrstr,istrip,nchp,  
      1                  qdiag(1,1,idlwdtc,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(idhdtc.gt.0) then  
       call paste2grd(dhsdtc,igrd,chfrstr,istrip,nchp,  
      1                  qdiag(1,1,idhdtc,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(idedtc.gt.0) then  
       call paste2grd(dedtc,igrd,chfrstr,istrip,nchp,  
      1                 qdiag(1,1,idedtc,bi,bj),ijall,1,nn,.false.)  
745        endif        endif
746        if(idhdqa.gt.0) then        if(diagnostics_is_on('GREEN   ',myid) ) then
747        call paste2grd(dhsdqa,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (grnstrip,igrd,chfrstr,istrip,nchp,nn,
748       1                  qdiag(1,1,idhdqa,bi,bj),ijall,1,nn,.false.)       . .false., 'GREEN   ', 1, 1, bi, bj, myid)
749          endif
750          if(diagnostics_is_on('LAI     ',myid) ) then
751          call diag_vegtile_fill (laistrip,igrd,chfrstr,istrip,nchp,nn,
752         . .false., 'LAI     ', 1, 1, bi, bj, myid)
753          endif
754          if(diagnostics_is_on('PARDR   ',myid) ) then
755          call diag_vegtile_fill (pardr,igrd,chfrstr,istrip,nchp,nn,
756         . .false., 'PARDR   ', 1, 1, bi, bj, myid)
757          endif
758          if(diagnostics_is_on('PARDF   ',myid) ) then
759          call diag_vegtile_fill (pardf,igrd,chfrstr,istrip,nchp,nn,
760         . .false., 'PARDF   ', 1, 1, bi, bj, myid)
761          endif
762          if(diagnostics_is_on('DLWDTC  ',myid) ) then
763          call diag_vegtile_fill (blwrad,igrd,chfrstr,istrip,nchp,nn,
764         . .false., 'DLWDTC  ', 1, 1, bi, bj, myid)
765          endif
766          if(diagnostics_is_on('DHDTC   ',myid) ) then
767          call diag_vegtile_fill (dhsdtc,igrd,chfrstr,istrip,nchp,nn,
768         . .false., 'DHDTC   ', 1, 1, bi, bj, myid)
769          endif
770          if(diagnostics_is_on('DEDTC   ',myid) ) then
771          call diag_vegtile_fill (dedtc,igrd,chfrstr,istrip,nchp,nn,
772         . .false., 'DEDTC   ', 1, 1, bi, bj, myid)
773          endif
774          if(diagnostics_is_on('DHDQA   ',myid) ) then
775          call diag_vegtile_fill (dhsdqa,igrd,chfrstr,istrip,nchp,nn,
776         . .false., 'DHDQA   ', 1, 1, bi, bj, myid)
777          endif
778          if(diagnostics_is_on('DEDQA   ',myid) ) then
779          call diag_vegtile_fill (dedqa,igrd,chfrstr,istrip,nchp,nn,
780         . .false., 'DEDQA   ', 1, 1, bi, bj, myid)
781          endif
782          if(diagnostics_is_on('LWGDOWN ',myid) ) then
783          call diag_vegtile_fill (hlwdwn,igrd,chfrstr,istrip,nchp,nn,
784         . .false., 'LWGDOWN ', 1, 1, bi, bj, myid)
785        endif        endif
786        if(idedqa.gt.0) then  C**********************************************************************
       call paste2grd(dedqa,igrd,chfrstr,istrip,nchp,  
      1                 qdiag(1,1,idedqa,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(ilwgdown.gt.0) then  
       call paste2grd(hlwdwn,igrd,chfrstr,istrip,nchp,  
      1                  qdiag(1,1,ilwgdown,bi,bj),ijall,1,nn,.false.)  
       endif  
 c**********************************************************************  
787    
788        if(nland.gt.0)then        if(nland.gt.0)then
789    
790         call tile (         call tile (
791       I   nland, timstp, types, rainls, raincon,  newsnow,  wspeed,       I   nland, timstp, types, rainls, raincon,  newsnow,  wspeed,
792       I   eturb,  dedqa,  dedtc,  hsturb, dhsdqa, dhsdtc,       I   eturb,  dedqa,  dedtc,  hsturb, dhsdqa, dhsdtc,
# Line 778  c*************************************** Line 805  c***************************************
805        if( nice.gt.0 ) then        if( nice.gt.0 ) then
806         call seaice ( nocean, timstp,     hice,         call seaice ( nocean, timstp,     hice,
807       .               eturb(nland+1),    dedtc(nland+1),       .               eturb(nland+1),    dedtc(nland+1),
808       .              hsturb(nland+1),   dhsdtc(nland+1),       .              hsturb(nland+1),   dhsdtc(nland+1),
809       .             qicestr(nland+1), dqicestr(nland+1),       .             qicestr(nland+1), dqicestr(nland+1),
810       .               swnet(nland+1),  lwstrip(nland+1), blwrad(nland+1),       .               swnet(nland+1),  lwstrip(nland+1), blwrad(nland+1),
811       .          pke(nland+1,nlay+1),    icest(nland+1),       .          pke(nland+1,nlay+1),    icest(nland+1),
812       .                  tc(nland+1),       qa(nland+1) )       .                  tc(nland+1),       qa(nland+1) )
813        endif        endif
814    
815  c***********************************************************************  C***********************************************************************
816  c   diagnostics: fill arrays for lsm output fields  C   diagnostics: fill arrays for lsm output fields
817  c***********************************************************************  C***********************************************************************
818    
819        if(irunoff.gt.0) then        if(diagnostics_is_on('RUNOFF  ',myid) ) then
820        call paste2grd(runoff,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (runoff,igrd,chfrstr,istrip,nchp,nn,
821       1                  qdiag(1,1,irunoff,bi,bj),ijall,1,nn,.false.)       . .false., 'RUNOFF  ', 1, 1, bi, bj, myid)
822        endif        endif
823        if(ifwsoil.gt.0) then        if(diagnostics_is_on('FWSOIL  ',myid) ) then
824        call paste2grd(fwsoil,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (fwsoil,igrd,chfrstr,istrip,nchp,nn,
825       1                  qdiag(1,1,ifwsoil,bi,bj),ijall,1,nn,.false.)       . .false., 'FWSOIL  ', 1, 1, bi, bj, myid)
826        endif        endif
827        if(igdrain.gt.0) then        if(diagnostics_is_on('GDRAIN  ',myid) ) then
828        call paste2grd(gdrain,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (gdrain,igrd,chfrstr,istrip,nchp,nn,
829       1                  qdiag(1,1,igdrain,bi,bj),ijall,1,nn,.false.)       . .false., 'GDRAIN  ', 1, 1, bi, bj, myid)
830        endif        endif
831        if(isnowmelt.gt.0) then        if(diagnostics_is_on('SNOWMELT',myid) ) then
832        call paste2grd(smelt,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (smelt,igrd,chfrstr,istrip,nchp,nn,
833       1                 qdiag(1,1,isnowmelt,bi,bj),ijall,1,nn,.false.)       . .false., 'SNOWMELT', 1, 1, bi, bj, myid)
834        endif        endif
835        if(ieveg.gt.0) then        if(diagnostics_is_on('EVEG    ',myid) ) then
836        call paste2grd(eveg,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (eveg,igrd,chfrstr,istrip,nchp,nn,
837       1                    qdiag(1,1,ieveg,bi,bj),ijall,1,nn,.false.)       . .false., 'EVEG    ', 1, 1, bi, bj, myid)
838        endif        endif
839        if(iesnow.gt.0) then        if(diagnostics_is_on('ESNOW   ',myid) ) then
840        call paste2grd(esno,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (esno,igrd,chfrstr,istrip,nchp,nn,
841       1                    qdiag(1,1,iesnow,bi,bj),ijall,1,nn,.false.)       . .false., 'ESNOW   ', 1, 1, bi, bj, myid)
842        endif        endif
843        if(iesoil.gt.0) then        if(diagnostics_is_on('ESOIL   ',myid) ) then
844        call paste2grd(esoi,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (esoi,igrd,chfrstr,istrip,nchp,nn,
845       1                    qdiag(1,1,iesoil,bi,bj),ijall,1,nn,.false.)       . .false., 'ESOIL   ', 1, 1, bi, bj, myid)
846        endif        endif
847        if(ieresv.gt.0) then        if(diagnostics_is_on('ERESV   ',myid) ) then
848        call paste2grd(eint,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (eint,igrd,chfrstr,istrip,nchp,nn,
849       1                    qdiag(1,1,ieresv,bi,bj),ijall,1,nn,.false.)       . .false., 'ERESV   ', 1, 1, bi, bj, myid)
850        endif        endif
851        if(ievpot.gt.0) then        if(diagnostics_is_on('EVPOT   ',myid) ) then
852        call paste2grd(evpot,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (evpot,igrd,chfrstr,istrip,nchp,nn,
853       1                     qdiag(1,1,ievpot,bi,bj),ijall,1,nn,.false.)       . .false., 'EVPOT   ', 1, 1, bi, bj, myid)
854        endif        endif
855        if(idtc.gt.0) then        if(diagnostics_is_on('DTC     ',myid) ) then
856        call paste2grd(strdg1,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (strdg1,igrd,chfrstr,istrip,nchp,nn,
857       1                      qdiag(1,1,idtc,bi,bj),ijall,1,nn,.false.)       . .false., 'DTC     ', 1, 1, bi, bj, myid)
858        endif        endif
859        if(idqc.gt.0) then        if(diagnostics_is_on('DQC     ',myid) ) then
860        call paste2grd(strdg2,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (strdg2,igrd,chfrstr,istrip,nchp,nn,
861       1                      qdiag(1,1,idqc,bi,bj),ijall,1,nn,.false.)       . .false., 'DQC     ', 1, 1, bi, bj, myid)
862        endif        endif
863        if(itcdtc.gt.0) then        if(diagnostics_is_on('TCDTC   ',myid) ) then
864        call paste2grd(strdg3,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (strdg3,igrd,chfrstr,istrip,nchp,nn,
865       1                      qdiag(1,1,itcdtc,bi,bj),ijall,1,nn,.false.)       . .false., 'TCDTC   ', 1, 1, bi, bj, myid)
866        endif        endif
867        if(iraddtc.gt.0) then        if(diagnostics_is_on('RADDTC  ',myid) ) then
868        call paste2grd(strdg4,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (strdg4,igrd,chfrstr,istrip,nchp,nn,
869       1                      qdiag(1,1,iraddtc,bi,bj),ijall,1,nn,.false.)       . .false., 'RADDTC  ', 1, 1, bi, bj, myid)
870        endif        endif
871        if(isensdtc.gt.0) then        if(diagnostics_is_on('SENSDTC ',myid) ) then
872        call paste2grd(strdg5,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (strdg5,igrd,chfrstr,istrip,nchp,nn,
873       1                     qdiag(1,1,isensdtc,bi,bj),ijall,1,nn,.false.)       . .false., 'SENSDTC ', 1, 1, bi, bj, myid)
874        endif        endif
875        if(ilatdtc.gt.0) then        if(diagnostics_is_on('LATDTC  ',myid) ) then
876        call paste2grd(strdg6,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (strdg6,igrd,chfrstr,istrip,nchp,nn,
877       1                      qdiag(1,1,ilatdtc,bi,bj),ijall,1,nn,.false.)       . .false., 'LATDTC  ', 1, 1, bi, bj, myid)
878        endif        endif
879        if(itddtc.gt.0) then        if(diagnostics_is_on('TDDTC   ',myid) ) then
880        call paste2grd(strdg7,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (strdg7,igrd,chfrstr,istrip,nchp,nn,
881       1                      qdiag(1,1,itddtc,bi,bj),ijall,1,nn,.false.)       . .false., 'TDDTC   ', 1, 1, bi, bj, myid)
882        endif        endif
883        if(iqcdtc.gt.0) then        if(diagnostics_is_on('QCDTC   ',myid) ) then
884        call paste2grd(strdg8,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (strdg8,igrd,chfrstr,istrip,nchp,nn,
885       1                      qdiag(1,1,iqcdtc,bi,bj),ijall,1,nn,.false.)       . .false., 'QCDTC   ', 1, 1, bi, bj, myid)
886        endif        endif
887  c***********************************************************************  C***********************************************************************
888    
889        if( ndlsm.gt.1 ) then        if( ndlsm.gt.1 ) then
890        call pstbitint(types,qdiaglsm(1,1),istrip,nchp,nchptot,1,nn)        call pstbitint(types,qdiaglsm(1,1),istrip,nchp,nchptot,1,nn)
# Line 909  c     call pstbmpit(igrdstr,qdiaglsm(1,5 Line 936  c     call pstbmpit(igrdstr,qdiaglsm(1,5
936        call pastit (capacity,capac,istrip,nchp,nchptot,1,nn)        call pastit (capacity,capac,istrip,nchp,nchptot,1,nn)
937        call pastit (snowdepth,snodep,istrip,nchp,nchptot,1,nn)        call pastit (snowdepth,snodep,istrip,nchp,nchptot,1,nn)
938    
939  c**********************************************************************  C**********************************************************************
940  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
941  c**********************************************************************  C**********************************************************************
942    
943        do i =1,istrip        do i =1,istrip
944        th(i,nlay+1) = tc(i) / pke(i,nlay+1)        th(i,nlay+1) = tc(i) / pke(i,nlay+1)
# Line 937  c*************************************** Line 964  c***************************************
964         stqflux(i,L) = transsh(i,L) * (sh(i,L+1)-sh(i,L))         stqflux(i,L) = transsh(i,L) * (sh(i,L+1)-sh(i,L))
965        enddo        enddo
966        enddo        enddo
967    
968  c tendency updates  C tendency updates
969  c ----------------  C ----------------
970        do  l=1,nlay        do  l=1,nlay
971        call strip2tile(uz(1,1,l),igrd,tmpstrip,nchp,ijall,        call strip2tile(uz(1,1,l),igrd,tmpstrip,nchp,ijall,
972       1                                                 istrip,1,nn)       1                                                 istrip,1,nn)
# Line 960  c ---------------- Line 987  c ----------------
987        do i =1,istrip        do i =1,istrip
988        tends(i) = ( th(i,l)-tmpstrip(i) )        tends(i) = ( th(i,l)-tmpstrip(i) )
989        enddo        enddo
990    
991        call pastit (tends,dtmsc(1,l),istrip,nchp,nchptot,1,nn)        call pastit (tends,dtmsc(1,l),istrip,nchp,nchptot,1,nn)
992    
993        call strip2tile(qz(1,1,l,1),igrd,tmpstrip,nchp,ijall,        call strip2tile(qz(1,1,l,1),igrd,tmpstrip,nchp,ijall,
# Line 967  c ---------------- Line 995  c ----------------
995        do i =1,istrip        do i =1,istrip
996        tends(i) = ( sh(i,l)-tmpstrip(i) )        tends(i) = ( sh(i,l)-tmpstrip(i) )
997        enddo        enddo
998    
999        call pastit (tends,dqmsc(1,l,1),istrip,nchp,nchptot,1,nn)        call pastit (tends,dqmsc(1,l,1),istrip,nchp,nchptot,1,nn)
1000    
1001        do nt = 1,ntracers-ptracers  c     do nt = 1,ntracers-ptracers
1002         call strip2tile(qz(1,1,L,ptracers+nt),igrd,tmpstrip,nchp,  c      call strip2tile(qz(1,1,L,ptracers+nt),igrd,tmpstrip,nchp,
1003       1                                             ijall,istrip,1,nn)  c    1                                             ijall,istrip,1,nn)
1004         do i =1,istrip  c      do i =1,istrip
1005          tends(i) = ( tracers(i,L,nt)-tmpstrip(i) )  c       tends(i) = ( tracers(i,L,nt)-tmpstrip(i) )
1006         enddo  c      enddo
1007         call pastit(tends,dqmsc(1,L,ptracers+nt),istrip,nchp,  c      call pastit(tends,dqmsc(1,L,ptracers+nt),istrip,nchp,
1008       .                                     nchptot,1,nn)  c    .                                     nchptot,1,nn)
1009        enddo  c     enddo
1010    
1011        enddo        enddo
1012    C *********************************************************************
1013    C **** increment diagnostic arrays for quantities saved in trbflx
1014    C *********************************************************************
1015    
1016  c *********************************************************************  C note: the order, logic, and scaling of the heat and moisture flux
1017  c **** increment diagnostic arrays for quantities saved in trbflx  C       diagnostics is critical!
1018  c *********************************************************************  C ------------------------------
   
 c note: the order, logic, and scaling of the heat and moisture flux  
 c       diagnostics is critical!  
 c ------------------------------  
1019    
1020        if(ievap.gt.0) then        if(diagnostics_is_on('EVAP    ',myid) ) then
1021        do i=1,istrip        do i=1,istrip
1022        tmpstrip(i) = stqflux(i,nlay) * 86400        tmpstrip(i) = stqflux(i,nlay) * 86400
1023        enddo        enddo
1024        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
1025       1                    qdiag(1,1,ievap,bi,bj),ijall,1,nn,.false.)       . .false., 'EVAP    ', 1, 1, bi, bj, myid)
1026        endif        endif
1027        if(ieflux.gt.0) then        if(diagnostics_is_on('EFLUX   ',myid) ) then
1028        do i=1,istrip        do i=1,istrip
1029        tmpstrip(i) = stqflux(i,nlay) * alhl        tmpstrip(i) = stqflux(i,nlay) * alhl
1030        enddo        enddo
1031        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
1032       1                    qdiag(1,1,ieflux,bi,bj),ijall,1,nn,.false.)       . .false., 'EFLUX   ', 1, 1, bi, bj, myid)
1033        endif        endif
1034        if(ihflux.gt.0) then        if(diagnostics_is_on('HFLUX   ',myid) ) then
1035        do i=1,istrip        do i=1,istrip
1036        tmpstrip(i) = sttflux(i,nlay) * cp        tmpstrip(i) = sttflux(i,nlay) * cp
1037        enddo        enddo
1038        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
1039       1                    qdiag(1,1,ihflux,bi,bj),ijall,1,nn,.false.)       . .false., 'HFLUX   ', 1, 1, bi, bj, myid)
1040        endif        endif
1041        if(ituflux.gt.0) then        if(diagnostics_is_on('TUFLUX  ',myid) ) then
1042        call paste2grd(stuflux,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (stuflux,igrd,chfrstr,istrip,nchp,nn,
1043       1                   qdiag(1,1,ituflux,bi,bj),ijall,nlay,nn,.false.)       . .false., 'TUFLUX  ', 0, nlay, bi, bj, myid)
1044        endif        endif
1045        if(itvflux.gt.0) then        if(diagnostics_is_on('TVFLUX  ',myid) ) then
1046        call paste2grd(stvflux,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (stvflux,igrd,chfrstr,istrip,nchp,nn,
1047       1                   qdiag(1,1,itvflux,bi,bj),ijall,nlay,nn,.false.)       . .false., 'TVFLUX  ', 0, nlay, bi, bj, myid)
1048        endif        endif
1049        if(ittflux.gt.0) then        if(diagnostics_is_on('TTFLUX  ',myid) ) then
1050        do l=1,nlay        do l=1,nlay
1051        do i=1,istrip        do i=1,istrip
1052        sttflux(i,l) = sttflux(i,l) * cp        sttflux(i,l) = sttflux(i,l) * cp
1053        enddo        enddo
1054        enddo        enddo
1055        call paste2grd(sttflux,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (sttflux,igrd,chfrstr,istrip,nchp,nn,
1056       1                   qdiag(1,1,ittflux,bi,bj),ijall,nlay,nn,.false.)       . .false., 'TTFLUX  ', 0, nlay, bi, bj, myid)
1057        endif        endif
1058        if(itqflux.gt.0) then        if(diagnostics_is_on('TQFLUX  ',myid) ) then
1059        do l=1,nlay        do l=1,nlay
1060        do i=1,istrip        do i=1,istrip
1061        stqflux(i,l) = stqflux(i,l) * alhl        stqflux(i,l) = stqflux(i,l) * alhl
1062        enddo        enddo
1063        enddo        enddo
1064        call paste2grd(stqflux,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (stqflux,igrd,chfrstr,istrip,nchp,nn,
1065       1                   qdiag(1,1,itqflux,bi,bj),ijall,nlay,nn,.false.)       . .false., 'TQFLUX  ', 0, nlay, bi, bj, myid)
1066        endif        endif
1067        if(iri.gt.0) call paste2grd(sri,igrd,chfrstr,istrip,nchp,        if(diagnostics_is_on('RI      ',myid) ) then
1068       1                       qdiag(1,1,iri,bi,bj),ijall,nlay,nn,.false.)        call diag_vegtile_fill (sri,igrd,chfrstr,istrip,nchp,nn,
1069        if(ikh.gt.0) call paste2grd(skh,igrd,chfrstr,istrip,nchp,       . .false., 'RI      ', 0, nlay, bi, bj, myid)
1070       1                       qdiag(1,1,ikh,bi,bj),ijall,nlay,nn,.false.)        endif
1071        if(ikm.gt.0) call paste2grd(skm,igrd,chfrstr,istrip,nchp,        if(diagnostics_is_on('KH      ',myid) ) then
1072       1                       qdiag(1,1,ikm,bi,bj),ijall,nlay,nn,.false.)        call diag_vegtile_fill (skh,igrd,chfrstr,istrip,nchp,nn,
1073        if(ict.gt.0) then       . .false., 'KH      ', 0, nlay, bi, bj, myid)
1074        call paste2grd(sct,igrd,chfrstr,istrip,nchp,        endif
1075       1                   qdiag(1,1,ict,bi,bj),ijall,1,nn,.false.)        if(diagnostics_is_on('KM      ',myid) ) then
1076        endif        call diag_vegtile_fill (skm,igrd,chfrstr,istrip,nchp,nn,
1077        if(icu.gt.0) then       . .false., 'KM      ', 0, nlay, bi, bj, myid)
1078        call paste2grd(scu,igrd,chfrstr,istrip,nchp,        endif
1079       1                   qdiag(1,1,icu,bi,bj),ijall,1,nn,.false.)        if(diagnostics_is_on('CT      ',myid) ) then
1080        endif        call diag_vegtile_fill (sct,igrd,chfrstr,istrip,nchp,nn,
1081        if(iwinds.gt.0) then       . .false., 'CT      ', 1, 1, bi, bj, myid)
1082        call paste2grd(swinds,igrd,chfrstr,istrip,nchp,        endif
1083       1                      qdiag(1,1,iwinds,bi,bj),ijall,1,nn,.false.)        if(diagnostics_is_on('CU      ',myid) ) then
1084        endif        call diag_vegtile_fill (scu,igrd,chfrstr,istrip,nchp,nn,
1085        if(iuflux.gt.0) then       . .false., 'CU      ', 1, 1, bi, bj, myid)
1086        call paste2grd(stuflux(1,nlay),igrd,chfrstr,istrip,nchp,        endif
1087       1                       qdiag(1,1,iuflux,bi,bj),ijall,1,nn,.false.)        if(diagnostics_is_on('WINDS   ',myid) ) then
1088        endif        call diag_vegtile_fill (swinds,igrd,chfrstr,istrip,nchp,nn,
1089        if(ivflux.gt.0) then       . .false., 'WINDS   ', 1, 1, bi, bj, myid)
1090        call paste2grd(stvflux(1,nlay),igrd,chfrstr,istrip,nchp,        endif
1091       1                       qdiag(1,1,ivflux,bi,bj),ijall,1,nn,.false.)        if(diagnostics_is_on('UFLUX   ',myid) ) then
1092        endif        call diag_vegtile_fill (stuflux(1,nlay),igrd,chfrstr,istrip,nchp,
1093        if(iustar.gt.0) then       . nn,.false., 'UFLUX   ', 1, 1, bi, bj, myid)
1094        call paste2grd(sustar,igrd,chfrstr,istrip,nchp,        endif
1095       1                      qdiag(1,1,iustar,bi,bj),ijall,1,nn,.false.)        if(diagnostics_is_on('VFLUX   ',myid) ) then
1096        endif        call diag_vegtile_fill (stvflux(1,nlay),igrd,chfrstr,istrip,nchp,
1097        if(iz0.gt.0) then       . nn,.false., 'VFLUX   ', 1, 1, bi, bj, myid)
1098        call paste2grd(sz0,igrd,chfrstr,istrip,nchp,        endif
1099       1                   qdiag(1,1,iz0,bi,bj),ijall,1,nn,.false.)        if(diagnostics_is_on('USTAR   ',myid) ) then
1100        endif        call diag_vegtile_fill (sustar,igrd,chfrstr,istrip,nchp,nn,
1101        if(ifrqtrb.gt.0) then       . .false., 'USTAR   ', 1, 1, bi, bj, myid)
1102        call paste2grd(frqtrb,igrd,chfrstr,istrip,nchp,        endif
1103       1                      qdiag(1,1,ifrqtrb,bi,bj),ijall,1,nn,.false.)        if(diagnostics_is_on('Z0      ',myid) ) then
1104        endif        call diag_vegtile_fill (sz0,igrd,chfrstr,istrip,nchp,nn,
1105        if(ipbl.gt.0) then       . .false., 'Z0      ', 1, 1, bi, bj, myid)
1106        call paste2grd(pbldpth,igrd,chfrstr,istrip,nchp,        endif
1107       1                       qdiag(1,1,ipbl,bi,bj),ijall,1,nn,.false.)        if(diagnostics_is_on('FRQTRB  ',myid) ) then
1108        endif        call diag_vegtile_fill (frqtrb,igrd,chfrstr,istrip,nchp,nn,
1109        if(iu2m.gt.0) then       . .false., 'FRQTRB  ', 0, nlay-1, bi, bj, myid)
1110        call paste2grd(stu2m,igrd,chfrstr,istrip,nchp,        endif
1111       1                     qdiag(1,1,iu2m,bi,bj),ijall,1,nn,.true.)        if(diagnostics_is_on('PBL     ',myid) ) then
1112        endif        call diag_vegtile_fill (pbldpth,igrd,chfrstr,istrip,nchp,nn,
1113        if(iv2m.gt.0) then       . .false., 'PBL     ', 1, 1, bi, bj, myid)
1114        call paste2grd(stv2m,igrd,chfrstr,istrip,nchp,        endif
1115       1                     qdiag(1,1,iv2m,bi,bj),ijall,1,nn,.true.)        if(diagnostics_is_on('U2M     ',myid) ) then
1116        endif        call diag_vegtile_fill (stu2m,igrd,chfrstr,istrip,nchp,nn,
1117        if(it2m.gt.0) then       . .false., 'U2M     ', 1, 1, bi, bj, myid)
1118        call paste2grd(stt2m,igrd,chfrstr,istrip,nchp,        endif
1119       1                     qdiag(1,1,it2m,bi,bj),ijall,1,nn,.true.)        if(diagnostics_is_on('V2M     ',myid) ) then
1120          call diag_vegtile_fill (stv2m,igrd,chfrstr,istrip,nchp,nn,
1121         . .false., 'V2M     ', 1, 1, bi, bj, myid)
1122          endif
1123          if(diagnostics_is_on('T2M     ',myid) ) then
1124          call diag_vegtile_fill (stt2m,igrd,chfrstr,istrip,nchp,nn,
1125         . .false., 'T2M     ', 1, 1, bi, bj, myid)
1126        endif        endif
1127        if(iq2m.gt.0) then        if(diagnostics_is_on('Q2M     ',myid) ) then
1128        do i=1,istrip        do i=1,istrip
1129           if( stq2m(i).ne.undef ) then           if( stq2m(i).ne.undef ) then
1130               tmpstrip(i) = stq2m(i) * 1000               tmpstrip(i) = stq2m(i) * 1000
# Line 1098  c ------------------------------ Line 1132  c ------------------------------
1132               tmpstrip(i) = undef               tmpstrip(i) = undef
1133           endif           endif
1134        enddo        enddo
1135        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
1136       1                     qdiag(1,1,iq2m,bi,bj),ijall,1,nn,.true.)       . .false., 'Q2M     ', 1, 1, bi, bj, myid)
1137        endif        endif
1138        if(iu10m.gt.0) then        if(diagnostics_is_on('U10M    ',myid) ) then
1139        call paste2grd(stu10m,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (stu10m,igrd,chfrstr,istrip,nchp,nn,
1140       1                      qdiag(1,1,iu10m,bi,bj),ijall,1,nn,.false.)       . .false., 'U10M    ', 1, 1, bi, bj, myid)
1141        endif        endif
1142        if(iv10m.gt.0) then        if(diagnostics_is_on('V10M    ',myid) ) then
1143        call paste2grd(stv10m,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (stv10m,igrd,chfrstr,istrip,nchp,nn,
1144       1                      qdiag(1,1,iv10m,bi,bj),ijall,1,nn,.false.)       . .false., 'V10M    ', 1, 1, bi, bj, myid)
1145        endif        endif
1146        if(it10m.gt.0) then        if(diagnostics_is_on('T10M    ',myid) ) then
1147        call paste2grd(stt10m,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (stt10m,igrd,chfrstr,istrip,nchp,nn,
1148       1                      qdiag(1,1,it10m,bi,bj),ijall,1,nn,.false.)       . .false., 'T10M    ', 1, 1, bi, bj, myid)
1149        endif        endif
1150        if(iq10m.gt.0) then        if(diagnostics_is_on('Q10M    ',myid) ) then
1151        do i=1,istrip        do i=1,istrip
1152           if( stq10m(i).ne.undef ) then           if( stq10m(i).ne.undef ) then
1153               tmpstrip(i) = stq10m(i) * 1000               tmpstrip(i) = stq10m(i) * 1000
# Line 1121  c ------------------------------ Line 1155  c ------------------------------
1155               tmpstrip(i) = undef               tmpstrip(i) = undef
1156           endif           endif
1157        enddo        enddo
1158        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
1159       1                      qdiag(1,1,iq10m,bi,bj),ijall,1,nn,.false.)       . .false., 'Q10M    ', 1, 1, bi, bj, myid)
1160        endif        endif
1161    
1162  c**********************************************************************  C**********************************************************************
1163  c   more diagnostics: land surface model parameters  C   more diagnostics: land surface model parameters
1164  c**********************************************************************  C**********************************************************************
   
       if(itdeep.gt.0)call paste2grd(td,igrd,chfrstr,istrip,nchp,  
      .                     qdiag(1,1,itdeep,bi,bj),ijall,1,nn,.false.)  
       if(iqcanopy .gt.0)call paste2grd(qa,igrd,chfrstr,istrip,nchp,  
      .                  qdiag(1,1,iqcanopy,bi,bj) ,ijall,1,nn,.false.)  
       if(ismshal  .gt.0)call paste2grd(swet1,igrd,chfrstr,istrip,nchp,  
      .                  qdiag(1,1,ismshal,bi,bj)  ,ijall,1,nn,.false.)  
       if(ismroot  .gt.0)call paste2grd(swet2,igrd,chfrstr,istrip,nchp,  
      .                  qdiag(1,1,ismroot,bi,bj)  ,ijall,1,nn,.false.)  
       if(ismdeep  .gt.0)call paste2grd(swet3,igrd,chfrstr,istrip,nchp,  
      .                  qdiag(1,1,ismdeep,bi,bj)  ,ijall,1,nn,.false.)  
       if(icapacity.gt.0)call paste2grd(capacity,igrd,chfrstr,istrip,  
      .             nchp,qdiag(1,1,icapacity,bi,bj),ijall,1,nn,.false.)  
       if(isnow.gt.0)call paste2grd(snowdepth,igrd,chfrstr,istrip,nchp,  
      .                  qdiag(1,1,isnow,bi,bj)    ,ijall,1,nn,.false.)  
1165    
1166  c**********************************************************************        if(diagnostics_is_on('TDEEP   ',myid) ) then
1167  c end regions loop        call diag_vegtile_fill (td,igrd,chfrstr,istrip,nchp,nn,
1168         . .false., 'TDEEP   ', 1, 1, bi, bj, myid)
1169          endif
1170          if(diagnostics_is_on('QCANOPY ',myid) ) then
1171          call diag_vegtile_fill (qa,igrd,chfrstr,istrip,nchp,nn,
1172         . .false., 'QCANOPY ', 1, 1, bi, bj, myid)
1173          endif
1174          if(diagnostics_is_on('SMSHAL  ',myid) ) then
1175          call diag_vegtile_fill (swet1,igrd,chfrstr,istrip,nchp,nn,
1176         . .false., 'SMSHAL  ', 1, 1, bi, bj, myid)
1177          endif
1178          if(diagnostics_is_on('SMROOT  ',myid) ) then
1179          call diag_vegtile_fill (swet2,igrd,chfrstr,istrip,nchp,nn,
1180         . .false., 'SMROOT  ', 1, 1, bi, bj, myid)
1181          endif
1182          if(diagnostics_is_on('SMDEEP  ',myid) ) then
1183          call diag_vegtile_fill (swet3,igrd,chfrstr,istrip,nchp,nn,
1184         . .false., 'SMDEEP  ', 1, 1, bi, bj, myid)
1185          endif
1186          if(diagnostics_is_on('CAPACITY',myid) ) then
1187          call diag_vegtile_fill (capacity,igrd,chfrstr,istrip,nchp,nn,
1188         . .false., 'CAPACITY', 1, 1, bi, bj, myid)
1189          endif
1190          if(diagnostics_is_on('SNOW    ',myid) ) then
1191          call diag_vegtile_fill (snowdepth,igrd,chfrstr,istrip,nchp,nn,
1192         . .false., 'SNOW    ', 1, 1, bi, bj, myid)
1193          endif
1194    C**********************************************************************
1195    C end regions loop
1196    
1197   2000 continue   2000 continue
1198    
1199  c**********************************************************************  C**********************************************************************
1200    
1201  c increment the counter for the accumulated fcc and qliq arrays  C increment the counter for the accumulated fcc and qliq arrays
1202  c ---------------------------------------------------------------------  C ---------------------------------------------------------------------
1203        imstturblw = imstturblw + 1        imstturblw = imstturblw + 1
1204        imstturbsw = imstturbsw + 1        imstturbsw = imstturbsw + 1
1205    
1206  c prevent ice or snow from melting  C prevent ice or snow from melting
1207  c ---------------------------------------------------------------------  C ---------------------------------------------------------------------
1208        do i =1,nchptot        do i =1,nchptot
1209        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
1210        enddo        enddo
1211    
1212  c Update tcanopy and ecanopy from the land points of the  C Update tcanopy and ecanopy from the points of the
1213  c           tground and qground arrays  C           tground and qground arrays
1214  c ---------------------------------------------------------------------  C ---------------------------------------------------------------------
1215        do i =1,nchplnd        do i =1,nchptot
1216         tcanopy(i) = tground(i)         tcanopy(i) = tground(i)
1217         ecanopy(i) = qground(i)         ecanopy(i) = qground(i)
1218        enddo        enddo
1219    
1220  C Initialize Tendencies and Couplings  C Initialize Tendencies and Couplings
1221  c -----------------------------------  C -----------------------------------
1222        do L = 1,nlay        do L = 1,nlay
1223         do i = 1,ijall         do i = 1,ijall
1224             duturb(i,1,L) = 0.             duturb(i,1,L) = 0.
# Line 1189  c ----------------------------------- Line 1236  c -----------------------------------
1236        enddo        enddo
1237    
1238  C Return Tendencies and Couplings to Grid Space  C Return Tendencies and Couplings to Grid Space
1239  c ---------------------------------------------  C ---------------------------------------------
1240        do l = 1,nlay        do l = 1,nlay
1241        call msc2grd(igrd,chfr,dumsc(1,L),nchp,nchptot,fracland,        call msc2grd(igrd,chfr,dumsc(1,L),nchp,nchptot,fracland,
1242       .                                              duturb(1,1,L),im,jm)       .                                              duturb(1,1,L),im,jm)
# Line 1204  c -------------------------------------- Line 1251  c --------------------------------------
1251        call msc2grd(igrd,chfr,  tke(1,L),nchp,nchptot,fracland,        call msc2grd(igrd,chfr,  tke(1,L),nchp,nchptot,fracland,
1252       .                                              qqgrid(1,1,L),im,jm)       .                                              qqgrid(1,1,L),im,jm)
1253    
1254        call msc2grd(igrd,chfr, fccmsc(1,L),nchp,nchptot,fracland,        call msc2grd(igrd,chfr, fccmsc(1,L),nchp,nchptot,fracland,
1255       .                                              fcctmp(1,1,L),im,jm)       .                                              fcctmp(1,1,L),im,jm)
1256        call msc2grd(igrd,chfr,qliqmsc(1,L),nchp,nchptot,fracland,        call msc2grd(igrd,chfr,qliqmsc(1,L),nchp,nchptot,fracland,
1257       .                                             qliqtmp(1,1,L),im,jm)       .                                             qliqtmp(1,1,L),im,jm)
1258        enddo        enddo
1259    
1260  c Reduce clouds from conditionally unstable layer  C Reduce clouds from conditionally unstable layer
1261  c -----------------------------------------------  C -----------------------------------------------
1262        call ctei ( tz,qz,fcctmp,qliqtmp,plz,pkz,pkht,im*jm,nlay )        call ctei ( tz,qz,fcctmp,qliqtmp,plz,pkz,pkht,im*jm,nlay )
1263    
1264  C Bumb Total Cloud Liquid Water and Fraction by Instantanious Values  C Bump Total Cloud Liquid Water and Fraction by Instantaneous Values
1265  c ------------------------------------------------------------------  C ------------------------------------------------------------------
1266        do l = 1,nlay        do l = 1,nlay
1267         do j=1,jm  c      do j=1,jm
1268         do i=1,im  c      do i=1,im
1269           fccavesw(i,j,L) =  fccavesw(i,j,L) +  fcctmp(i,j,L)  c        fccavesw(i,j,L) =  fccavesw(i,j,L) +  fcctmp(i,j,L)
1270           fccavelw(i,j,L) =  fccavelw(i,j,L) +  fcctmp(i,j,L)  c        fccavelw(i,j,L) =  fccavelw(i,j,L) +  fcctmp(i,j,L)
1271          qliqavelw(i,j,L) = qliqavelw(i,j,L) + qliqtmp(i,j,L)  c       qliqavelw(i,j,L) = qliqavelw(i,j,L) + qliqtmp(i,j,L)
1272          qliqavesw(i,j,L) = qliqavesw(i,j,L) + qliqtmp(i,j,L)  c       qliqavesw(i,j,L) = qliqavesw(i,j,L) + qliqtmp(i,j,L)
1273         enddo  c      enddo
1274         enddo  c      enddo
1275           do i = 1,ijall
1276         if (itrbfcc.gt.0) then           fccavesw(i,1,L) =  fccavesw(i,1,L) +  fcctmp(i,1,L)
1277         do j=1,jm           fccavelw(i,1,L) =  fccavelw(i,1,L) +  fcctmp(i,1,L)
1278         do i=1,im          qliqavelw(i,1,L) = qliqavelw(i,1,L) + qliqtmp(i,1,L)
1279         qdiag(i,j,itrbfcc+L-1,bi,bj) =  qdiag(i,j,itrbfcc+L-1,bi,bj) +          qliqavesw(i,1,L) = qliqavesw(i,1,L) + qliqtmp(i,1,L)
      .                                              fcctmp(i,j,L)  
        enddo  
        enddo  
        endif  
   
        if (itrbqliq.gt.0) then  
        do j=1,jm  
        do i=1,im  
        qdiag(i,j,itrbqliq+L-1,bi,bj)=qdiag(i,j,itrbqliq+L-1,bi,bj)+  
      .                                             qliqtmp(i,j,L)*1.e6  
        enddo  
1280         enddo         enddo
        endif  
       enddo  
1281    
1282    c      if(diagnostics_is_on('TRBFCC  ',myid) ) then
1283    c       call diagnostics_fill(fcctmp(i,j,L),'TRBFCC  ',L,1,3,bi,bj,myid)
1284    c      endif
1285    c      if(diagnostics_is_on('TRBQLIQ ',myid) ) then
1286    c       do j = 1,jm
1287    c       do i = 1,im
1288    c        tmpdiag(i,j) = qliqtmp(i,j,L)*1.e6
1289    c       enddo
1290    c       enddo
1291    c       call diagnostics_fill(tmpdiag,'TRBQLIQ ',L,1,3,bi,bj,myid)
1292    c      endif
1293          enddo
1294          tmpFac = 1.e6
1295          CALL DIAGNOSTICS_FILL(fcctmp,'TRBFCC  ',0,nlay,3,bi,bj,myid)
1296          CALL DIAGNOSTICS_SCALE_FILL( qliqtmp,tmpFac,1,
1297         &                             'TRBQLIQ ',0,nlay,3,bi,bj,myid)
1298  C**********************************************************************  C**********************************************************************
1299  C And some other variables to be transformed back to grid space:  C And some other variables to be transformed back to grid space:
1300  C Ground Temperature, snow depth and shallow layer ground wetness  C Ground Temperature, snow depth and shallow layer ground wetness
1301        do j = 1,jm  c     do j = 1,jm
1302         do i = 1,im  c      do i = 1,im
1303           tgz(i,j) = 0.  c        tgz(i,j) = 0.
1304         enddo  c      enddo
1305    c     enddo
1306          do i = 1,ijall
1307             tgz(i,1) = 0.
1308        enddo        enddo
1309        call msc2grd(igrd,chfr,tground ,nchp,nchptot,fracland,tgz ,im,jm)        call msc2grd(igrd,chfr,tground ,nchp,nchptot,fracland,tgz ,im,jm)
1310    
1311  c *********************************************************************  C *********************************************************************
1312  c **** increment diagnostic array for ground and surface temperatures,  C **** increment diagnostic array for ground and surface temperatures,
1313  c ***       ground temp tendency, and ground humidity  C ***       ground temp tendency, and ground humidity
1314  c *********************************************************************  C *********************************************************************
   
       if(itground.gt.0) then  
       do i =1,ijall  
       qdiag(i,1,itground,bi,bj) = qdiag(i,1,itground,bi,bj) + tgz(i,1)  
       enddo  
       endif  
   
       if(itcanopy.gt.0) then  
       do i =1,ijall  
       qdiag(i,1,itcanopy,bi,bj) = qdiag(i,1,itcanopy,bi,bj) + tgz(i,1)  
       enddo  
       endif  
   
       if(its.gt.0) then  
       do i =1,ijall  
       tmpstrip(i) = tz(i,1,nlay) * pkht(i,1,nlay)  
       enddo  
       do i =1,ijall  
       qdiag(i,1,its,bi,bj) = qdiag(i,1,its,bi,bj) + tmpstrip(i)  
       enddo  
       endif  
   
       if(idtg.gt.0) then  
       do i =1,ijall  
       qdiag(i,1,idtg,bi,bj) = qdiag(i,1,idtg,bi,bj) + tgz(i,1)  
       enddo  
       endif  
1315    
1316  c *********************************************************************  c     if(diagnostics_is_on('TGROUND ',myid) ) then
1317  c ****           increment diagnostic arrays for tendencies        ****         call diagnostics_fill(tgz,'TGROUND ',0,1,3,bi,bj,myid)
1318  c *********************************************************************  c     endif
1319        do L = 1,nlay  c     if(diagnostics_is_on('TCANOPY ',myid) ) then
1320           call diagnostics_fill(tgz,'TCANOPY ',0,1,3,bi,bj,myid)
1321         if(iturbu.gt.0) then  c     endif
        do i =1,ijall  
         qdiag(i,1,iturbu+l-1,bi,bj) =  qdiag(i,1,iturbu+l-1,bi,bj)  
      .                      + duturb(i,1,l) * atimstp * secday  
        enddo  
        endif  
1322    
1323         if(iturbv.gt.0) then        if(diagnostics_is_on('TS      ',myid) ) then
1324         do i =1,ijall  c      do j =1,jm
1325          qdiag(i,1,iturbv+l-1,bi,bj) =  qdiag(i,1,iturbv+l-1,bi,bj)  c      do i =1,im
1326       .                      + dvturb(i,1,l) * atimstp * secday  c       tmpdiag(i,j) = tz(i,j,nlay) * pkht(i,j,nlay)
1327    c      enddo
1328           do i = 1,ijall
1329            tmpdiag(i,1) = tz(i,1,nlay) * pkht(i,1,nlay)
1330         enddo         enddo
1331         endif         call diagnostics_fill(tmpdiag,'TS      ',0,1,3,bi,bj,myid)
1332          endif
1333    
1334         if(iturbq.gt.0) then  c     if(diagnostics_is_on('DTG     ',myid) ) then
1335         do i =1,ijall         call diagnostics_fill(tgz,'DTG     ',0,1,3,bi,bj,myid)
1336          qdiag(i,1,iturbq+l-1,bi,bj) =  qdiag(i,1,iturbq+l-1,bi,bj)  c     endif
      .                      + dqturb(i,1,l,1) * atimstp * secday * 1000  
        enddo  
        endif  
1337    
1338         if(iturbt.gt.0) then  C *********************************************************************
1339         do i =1,ijall  C ****           increment diagnostic arrays for tendencies        ****
1340          qdiag(i,1,iturbt+l-1,bi,bj) =  qdiag(i,1,iturbt+l-1,bi,bj)  C *********************************************************************
1341       .                   + dtturb(i,1,l) * pkz(i,1,l)*atimstp*secday        tmpFac = atimstp * secday
1342         enddo        CALL DIAGNOSTICS_SCALE_FILL(dvturb,tmpFac,1,
1343         &                            'TURBV   ',0,nlay,3,bi,bj,myid)
1344          CALL DIAGNOSTICS_SCALE_FILL(duturb,tmpFac,1,
1345         &                            'TURBU   ',0,nlay,3,bi,bj,myid)
1346          tmpFac = atimstp * secday * 1000.
1347          CALL DIAGNOSTICS_SCALE_FILL(dqturb,tmpFac,1,
1348         &                            'TURBQ   ',0,nlay,3,bi,bj,myid)
1349    
1350    c     do L = 1,nlay
1351    
1352    c      if(diagnostics_is_on('TURBV   ',myid) ) then
1353    c       do j = 1,jm
1354    c       do i = 1,im
1355    c        tmpdiag(i,j) = dvturb(i,j,l) * atimstp * secday
1356    c       enddo
1357    c       enddo
1358    c       call diagnostics_fill(tmpdiag,'TURBV   ',L,1,3,bi,bj,myid)
1359    c      endif
1360    
1361    c      if(diagnostics_is_on('TURBU   ',myid) ) then
1362    c       do j = 1,jm
1363    c       do i = 1,im
1364    c        tmpdiag(i,j) = duturb(i,j,l) * atimstp * secday
1365    c       enddo
1366    c       enddo
1367    c       call diagnostics_fill(tmpdiag,'TURBU   ',L,1,3,bi,bj,myid)
1368    c      endif
1369    
1370    c      if(diagnostics_is_on('TURBQ   ',myid) ) then
1371    c       do j = 1,jm
1372    c       do i = 1,im
1373    c        tmpdiag(i,j) = dqturb(i,j,l,1) * atimstp * secday * 1000.
1374    c       enddo
1375    c       enddo
1376    c       call diagnostics_fill(tmpdiag,'TURBQ   ',L,1,3,bi,bj,myid)
1377    c      endif
1378    
1379           if(diagnostics_is_on('TURBT   ',myid) ) then
1380    c       do j = 1,jm
1381    c       do i = 1,im
1382    c        tmpdiag(i,j) = dtturb(i,j,l) * pkz(i,j,l)*atimstp*secday
1383    c       enddo
1384    c       enddo
1385            do L = 1,nlay
1386             do i = 1,ijall
1387              tmpdiag(i,1) = dtturb(i,1,l) * pkz(i,1,l)*atimstp*secday
1388             enddo
1389             call diagnostics_fill(tmpdiag,'TURBT   ',L,1,3,bi,bj,myid)
1390            enddo
1391         endif         endif
1392    
1393        enddo  c     enddo
1394    
1395  c pi-weight the theta and moisture tendencies  C pi-weight the theta and moisture tendencies
1396  c -------------------------------------------  C -------------------------------------------
1397        do i =1,ijall        do i =1,ijall
1398        thtgz(i) = pz(i,1) * atimstp        thtgz(i) = pz(i,1) * atimstp
1399        enddo        enddo
# Line 1340  c -------------------------------------- Line 1410  c --------------------------------------
1410         enddo         enddo
1411        enddo        enddo
1412    
1413  c *********************************************************************  C *********************************************************************
1414  c ****   zero out the accumulating rainfall and snowfall arrays     ***  C ****   zero out the accumulating rainfall and snowfall arrays     ***
1415  c *********************************************************************  C *********************************************************************
1416    
1417        if( time_left.lt.timstp ) then        if( time_left.lt.timstp ) then
1418        do j = 1,jm  c     do j = 1,jm
1419        do i = 1,im  c     do i = 1,im
1420           rainlsp(i,j) = 0.  c        rainlsp(i,j) = 0.
1421          rainconv(i,j) = 0.  c       rainconv(i,j) = 0.
1422          snowfall(i,j) = 0.  c       snowfall(i,j) = 0.
1423        enddo  c     enddo
1424        enddo  c     enddo
1425        endif         do i =1,ijall
1426             rainlsp(i,1) = 0.
1427  c *********************************************************************          rainconv(i,1) = 0.
1428  c ****                bump diagnostic counters                      ***          snowfall(i,1) = 0.
1429  c *********************************************************************         enddo
1430          endif
       nturbt    = nturbt    + 1  
       nturbq    = nturbq    + 1  
       nturbu    = nturbu    + 1  
       nturbv    = nturbv    + 1  
       ntuflux   = ntuflux   + 1  
       ntvflux   = ntvflux   + 1  
       nttflux   = nttflux   + 1  
       ntqflux   = ntqflux   + 1  
       nwinds    = nwinds    + 1  
       nkm       = nkm       + 1  
       nkh       = nkh       + 1  
       nri       = nri       + 1  
       nct       = nct       + 1  
       ncu       = ncu       + 1  
       ntground  = ntground  + 1  
       nts       = nts       + 1  
       ndtg      = ndtg      + 1  
       nqg       = nqg       + 1  
       nqs       = nqs       + 1  
       nhflux    = nhflux    + 1  
       neflux    = neflux    + 1  
       nevap     = nevap     + 1  
       nuflux    = nuflux    + 1  
       nvflux    = nvflux    + 1  
       ndtsrf    = ndtsrf    + 1  
       nustar    = nustar    + 1  
       nz0       = nz0       + 1  
       nfrqtrb   = nfrqtrb   + 1  
       npbl      = npbl      + 1  
       nu2m      = nu2m      + 1  
       nv2m      = nv2m      + 1  
       nt2m      = nt2m      + 1  
       nq2m      = nq2m      + 1  
       nu10m     = nu10m     + 1  
       nv10m     = nv10m     + 1  
       nt10m     = nt10m     + 1  
       nq10m     = nq10m     + 1  
       ntcanopy  = ntcanopy  + 1  
       ntdeep    = ntdeep    + 1  
       nqcanopy  = nqcanopy  + 1  
       nsmshal   = nsmshal   + 1  
       nsmroot   = nsmroot   + 1  
       nsmdeep   = nsmdeep   + 1  
       nsnow     = nsnow     + 1  
       ncapacity = ncapacity + 1  
       nraincon  = nraincon  + 1  
       nrainlsp  = nrainlsp  + 1  
       nsnowfall = nsnowfall + 1  
       nrunoff   = nrunoff   + 1  
       nfwsoil   = nfwsoil   + 1  
       ngdrain   = ngdrain   + 1  
       nsnowmelt = nsnowmelt + 1  
       neresv    = neresv    + 1  
       nesoil    = nesoil    + 1  
       neveg     = neveg     + 1  
       nesnow    = nesnow    + 1  
       npardf    = npardf    + 1  
       npardr    = npardr    + 1  
       nlai      = nlai      + 1  
       ngreen    = ngreen    + 1  
       ndlwdtc   = ndlwdtc   + 1  
       ndhdtc    = ndhdtc    + 1  
       ndedtc    = ndedtc    + 1  
       nevpot    = nevpot    + 1  
       nlwgdown  = nlwgdown  + 1  
       ndhdqa    = ndhdqa    + 1  
       ndedqa    = ndedqa    + 1  
       ndtc      = ndtc      + 1  
       ndqc      = ndqc      + 1  
       ntcdtc    = ntcdtc    + 1  
       nraddtc   = nraddtc   + 1  
       nsensdtc  = nsensdtc  + 1  
       nlatdtc   = nlatdtc   + 1  
       ntddtc    = ntddtc    + 1  
       nqcdtc    = nqcdtc    + 1  
       ntrbqliq  = ntrbqliq  + 1  
       ntrbfcc   = ntrbfcc   + 1  
1431    
1432        return        return
1433        end        end
# Line 1442  c ************************************** Line 1435  c **************************************
1435       1 IWATER,Z0,tracers,ntrace,ntracedim,DTAU,ITRTRB,KMBG,KHBG,QBEG,       1 IWATER,Z0,tracers,ntrace,ntracedim,DTAU,ITRTRB,KMBG,KHBG,QBEG,
1436       2 TPROF,WU,WV,SRI,ET,EU,SWINDS,sustar,sz0,freqdg,pbldpth,       2 TPROF,WU,WV,SRI,ET,EU,SWINDS,sustar,sz0,freqdg,pbldpth,
1437       3 sct,scu,stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,stq10m,       3 sct,scu,stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,stq10m,
1438       4 irun,nlev,NYMD,NHMS,grav,cp,rgas,faceps,virtcon,undef,       4 irun,nlev,nltop,NYMD,NHMS,grav,cp,rgas,faceps,virtcon,undef,
1439       5 dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,       5 dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,
1440       6 hsturb,dhsdqa,dhsdtc,transth,transsh,       6 hsturb,dhsdqa,dhsdtc,transth,transsh,
1441       7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)       7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)
# Line 1468  C    PLKE          -         EDGE LEVEL Line 1461  C    PLKE          -         EDGE LEVEL
1461  C    DPSTR         -         PRESSURE INTERVALS  C    DPSTR         -         PRESSURE INTERVALS
1462  C    WATER         -         BIT ARRAY - '1' OVER OCEANS  C    WATER         -         BIT ARRAY - '1' OVER OCEANS
1463  C    Z0            -         SURFACE ROUGHNESS  C    Z0            -         SURFACE ROUGHNESS
1464  C    tracers       -         array of passive tracers  C    tracers       -         array of passive tracers
1465  C    ntrace        -         number of tracers to be diffused  C    ntrace        -         number of tracers to be diffused
1466  C    ntracedim     -         outer dimension of tracers array  C    ntracedim     -         outer dimension of tracers array
1467  C    DTAU          -         TIME CHANGE PER ITERATION OF TRBLFX  C    DTAU          -         TIME CHANGE PER ITERATION OF TRBFLX
1468  C    ITRTRB        -         NUMBER OF ITERATIONS OF TRBLFX  C    ITRTRB        -         NUMBER OF ITERATIONS OF TRBFLX
1469  C    KMBG          -         BACKGROUND VALUE OF MOMENTUM TRANSFER COEF  C    KMBG          -         BACKGROUND VALUE OF MOMENTUM TRANSFER COEF
1470  C    KHBG          -         BACKGROUND VALUE OF HEAT TRANSFER COEF  C    KHBG          -         BACKGROUND VALUE OF HEAT TRANSFER COEF
1471  C    NLEV          -         NUMBER OF ATMOSPHERIC LEVELS TO CALCULATE  C    NLEV          -         NUMBER OF ATMOSPHERIC LEVELS TO CALCULATE
1472    C    nltop         -         Top level allowed for turbulence
1473  C    QBEG          -         LOGICAL .TRUE. FOR INITIAL START OF GCM  C    QBEG          -         LOGICAL .TRUE. FOR INITIAL START OF GCM
1474  C    TPROF         -         LOGICAL .TRUE. TO CALCULATE PT BY PT DIAGS  C    TPROF         -         LOGICAL .TRUE. TO CALCULATE PT BY PT DIAGS
1475  C  C
# Line 1487  C*************************************** Line 1481  C***************************************
1481        implicit none        implicit none
1482    
1483  C Argument list declarations  C Argument list declarations
1484        integer nn,irun,nlev,ntrace,ntracedim,itrtrb,nhms,nymd        integer nn,irun,nlev,nltop,ntrace,ntracedim,itrtrb,nhms,nymd
1485        real TH(irun,NLEV+1),THV(irun,NLEV+1),SH(irun,NLEV+1)        _RL TH(irun,NLEV+1),THV(irun,NLEV+1),SH(irun,NLEV+1)
1486        real U(irun,NLEV+1),V(irun,NLEV+1),QQ(irun,NLEV)        _RL U(irun,NLEV+1),V(irun,NLEV+1),QQ(irun,NLEV)
1487        real PL(irun,NLEV),PLE(irun,NLEV+1),PLK(irun,NLEV)        _RL PL(irun,NLEV),PLE(irun,NLEV+1),PLK(irun,NLEV)
1488        real PLKE(irun,NLEV+1),DPSTR(irun,NLEV)        _RL PLKE(irun,NLEV+1),DPSTR(irun,NLEV)
1489        integer IWATER(irun)        integer IWATER(irun)
1490        real Z0(irun)        _RL Z0(irun)
1491        real tracers(irun,nlev+1,ntracedim)        _RL tracers(irun,nlev+1,ntracedim)
1492        real dtau,KMBG,KHBG        _RL dtau,KMBG,KHBG
1493        LOGICAL QBEG,TPROF        LOGICAL QBEG,TPROF
1494        real SWINDS(irun)        _RL SWINDS(irun)
1495        real SRI(irun,nlev), ET(irun,nlev)        _RL SRI(irun,nlev), ET(irun,nlev)
1496        real EU (irun,nlev)        _RL EU (irun,nlev)
1497        real WU(irun,nlev)        _RL WU(irun,nlev)
1498        real WV (irun,nlev), pbldpth(irun)        _RL WV (irun,nlev), pbldpth(irun)
1499        real sustar(irun), sz0(irun)        _RL sustar(irun), sz0(irun)
1500        real freqdg(irun,nlev-1)        _RL freqdg(irun,nlev-1)
1501        real sct(irun), scu(irun)        _RL sct(irun), scu(irun)
1502        real stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)        _RL stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)
1503        real stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)        _RL stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)
1504        real grav,cp,rgas,faceps,virtcon,undef        _RL grav,cp,rgas,faceps,virtcon,undef
1505        real eturb(irun),dedqa(irun),dedtc(irun)        _RL eturb(irun),dedqa(irun),dedtc(irun)
1506        real hsturb(irun),dhsdqa(irun),dhsdtc(irun)        _RL hsturb(irun),dhsdqa(irun),dhsdtc(irun)
1507        real dshdthg(irun,nlev),dthdthg(irun,nlev)        _RL dshdthg(irun,nlev),dthdthg(irun,nlev)
1508        real dshdshg(irun,nlev),dthdshg(irun,nlev)        _RL dshdshg(irun,nlev),dthdshg(irun,nlev)
1509        real transth(irun,nlev),transsh(irun,nlev)        _RL transth(irun,nlev),transsh(irun,nlev)
1510        real ctsave(irun),xxsave(irun),yysave(irun)        _RL ctsave(irun),xxsave(irun),yysave(irun)
1511        real zetasave(irun),xlsave(irun,nlev),khsave(irun,nlev)        _RL zetasave(irun),xlsave(irun,nlev),khsave(irun,nlev)
1512        real qliq(irun,nlev),turbfcc(irun,nlev)        _RL qliq(irun,nlev),turbfcc(irun,nlev)
1513    
1514  C Local Variables  C Local Variables
1515        real b1,b3,alpha,halpha,qqmin,qbustr        _RL b1,b3,alpha,halpha,qqmin,qbustr
1516        PARAMETER ( B1      =   16.6    )          PARAMETER ( B1      =   16.6    )
1517        PARAMETER ( B3      = 1. / B1  )          PARAMETER ( B3      = 1. / B1  )
1518        PARAMETER ( ALPHA   = 0.1       )        PARAMETER ( ALPHA   = 0.1       )
1519        PARAMETER ( HALPHA = ALPHA * 0.5 )        PARAMETER ( HALPHA = ALPHA * 0.5 )
1520        PARAMETER ( QQMIN  = 0.005      )        PARAMETER ( QQMIN  = 0.005      )
1521        PARAMETER ( QBUSTR = 2.550952   )        PARAMETER ( QBUSTR = 2.550952   )
1522        real    argmax, onethrd, z1pem25, b2, two        _RL    argmax, onethrd, z1pem25, b2, two
1523        PARAMETER (ARGMAX = 30.)        PARAMETER (ARGMAX = 30.)
1524        PARAMETER (ONETHRD = 1./3. )        PARAMETER (ONETHRD = 1./3. )
1525        PARAMETER (Z1PEM25 = 1.E-25)        PARAMETER (Z1PEM25 = 1.E-25)
1526        PARAMETER ( B2    =  10.1 )        PARAMETER ( B2    =  10.1 )
1527        PARAMETER ( two   =   2.0 )        PARAMETER ( two   =   2.0 )
1528    
1529        real AHS (irun), HS(irun)        _RL AHS (irun), HS(irun)
1530        real XX  (irun), YY(irun), CU(irun)        _RL XX  (irun), YY(irun), CU(irun)
1531        real CT(irun),  USTAR(irun)        _RL CT(irun),  USTAR(irun)
1532        real RIB(irun),   ZETA(irun), WS(irun)        _RL RIB(irun),   ZETA(irun), WS(irun)
1533        real DTHS(irun), DELTHS(irun)        _RL DTHS(irun), DELTHS(irun)
1534        real DTHL(irun), DELTHL(irun)        _RL DTHL(irun), DELTHL(irun)
1535        real RIBIN(irun),CUIN(irun)        _RL RIBIN(irun),CUIN(irun)
1536        real CTIN(irun),ZETAIN(irun)        _RL CTIN(irun),ZETAIN(irun)
1537        real USTARIN(irun),RHOSIN(irun),Z0IN(irun)        _RL USTARIN(irun),RHOSIN(irun),Z0IN(irun)
1538        real qqcolmin(irun),qqcolmax(irun),levpbl(irun)        _RL qqcolmin(irun),qqcolmax(irun),levpbl(irun)
1539    
1540        real ADZ1(irun,nlev), DZ1TMP(irun,nlev)        _RL ADZ1(irun,nlev), DZ1TMP(irun,nlev)
1541        real DZ3(irun,nlev), TEMP(irun,nlev)        _RL DZ3(irun,nlev), TEMP(irun,nlev)
1542        real DV(irun,nlev), DTHV(irun,nlev)        _RL DV(irun,nlev), DTHV(irun,nlev)
1543        real DPK(irun,nlev), STRT(irun,nlev)        _RL DPK(irun,nlev), STRT(irun,nlev)
1544        real DW2(irun,nlev), RI(irun,nlev)        _RL DW2(irun,nlev), RI(irun,nlev)
1545        real RHOZPK(irun,nlev), Q(irun,nlev)        _RL RHOZPK(irun,nlev), Q(irun,nlev)
1546        real RIINIT(irun,nlev), DU(irun,nlev)        _RL RIINIT(irun,nlev), DU(irun,nlev)
1547        real QQINIT(irun,nlev), RHOKDZ(irun,nlev)        _RL QQINIT(irun,nlev), RHOKDZ(irun,nlev)
1548        real RHODZ2(irun,nlev)        _RL RHODZ2(irun,nlev)
1549        REAL KM(irun,nlev), KH(irun,nlev)        _RL KM(irun,nlev), KH(irun,nlev)
1550    
1551        real DELTH  (irun,nlev+1), DELSH (irun,nlev+1)        _RL DELTH  (irun,nlev+1), DELSH (irun,nlev+1)
1552        real FLXFAC (irun,nlev+1)        _RL FLXFAC (irun,nlev+1)
1553        real FLXFPK (irun,nlev+1)        _RL FLXFPK (irun,nlev+1)
1554    
1555        real ADZ2   (irun,nlev-1), RHODZ1(irun,nlev-1)        _RL ADZ2   (irun,nlev-1), RHODZ1(irun,nlev-1)
1556        real VKZE   (irun,nlev-1), VKZM  (irun,nlev-1)        _RL VKZE   (irun,nlev-1), VKZM  (irun,nlev-1)
1557        real XL     (irun,nlev-1), QXLM  (irun,nlev-1)        _RL XL     (irun,nlev-1), QXLM  (irun,nlev-1)
1558        real QQE    (irun,nlev-1), QE    (irun,nlev-1)        _RL QQE    (irun,nlev-1), QE    (irun,nlev-1)
1559        real P3     (irun,nlev-1), XQ    (irun,nlev-1)        _RL P3     (irun,nlev-1), XQ    (irun,nlev-1)
1560        real XLDIAG (irun,nlev-1), FLXFCE(irun,nlev-1)        _RL XLDIAG (irun,nlev-1), FLXFCE(irun,nlev-1)
1561    
1562        LOGICAL FIRST,LAST        LOGICAL FIRST,LAST
1563        integer IBITSTB(irun,nlev),INTQ(irun,nlev)        integer IBITSTB(irun,nlev),INTQ(irun,nlev)
1564    
1565  C arrays for use by moist bouyancy calculation  C arrays for use by moist bouyancy calculation
1566  C -----------------  C -----------------
1567        real TL(irun,NLEV),DTH(irun,NLEV)        _RL TL(irun,NLEV),DTH(irun,NLEV)
1568        real DSH(irun,NLEV)        _RL DSH(irun,NLEV)
1569        real SHL(irun,NLEV)        _RL SHL(irun,NLEV)
1570        real AA(irun,NLEV),BB(irun,NLEV),SSDEV(irun,NLEV)        _RL AA(irun,NLEV),BB(irun,NLEV),SSDEV(irun,NLEV)
1571        real ARG(irun,NLEV),XXZETA(irun),QBYU(irun)        _RL ARG(irun,NLEV),XXZETA(irun),QBYU(irun)
1572        real SVAR(irun,NLEV),Q1M(irun,NLEV)        _RL SVAR(irun,NLEV),Q1M(irun,NLEV)
1573        real FCC(irun,NLEV)        _RL FCC(irun,NLEV)
1574        real BETAT(irun,NLEV),BETAW(irun,NLEV)        _RL BETAT(irun,NLEV),BETAW(irun,NLEV)
1575        real BETAL(irun,NLEV),BETAT1(irun,NLEV)        _RL BETAL(irun,NLEV),BETAT1(irun,NLEV)
1576        real BETAW1(irun,NLEV),SBAR(irun,NLEV)        _RL BETAW1(irun,NLEV),SBAR(irun,NLEV)
1577        real SHSAT(irun,NLEV)        _RL SHSAT(irun,NLEV)
1578    
1579  C Some space for variables to be used in called routines  C Some space for variables to be used in called routines
1580        logical LWATER        logical LWATER
1581        integer IVBITRIB(irun)        integer IVBITRIB(irun)
1582        real VHZ(irun)        _RL VHZ(irun)
1583        real VH0(irun)        _RL VH0(irun)
1584        real VPSIM(irun),VAPSIM(irun)        _RL VPSIM(irun),VAPSIM(irun)
1585        real VPSIG(irun),VPSIHG(irun)        _RL VPSIG(irun),VPSIHG(irun)
1586        real VTEMP(irun),VDZETA(irun)        _RL VTEMP(irun),VDZETA(irun)
1587        real VDZ0(irun),VDPSIM(irun)        _RL VDZ0(irun),VDPSIM(irun)
1588        real VDPSIH(irun),VZH(irun)        _RL VDPSIH(irun),VZH(irun)
1589        real VXX0(irun),VYY0(irun)        _RL VXX0(irun),VYY0(irun)
1590        real VAPSIHG(irun),VRIB1(irun),VWS1(irun)        _RL VAPSIHG(irun),VRIB1(irun),VWS1(irun)
1591        real VPSIH(irun),VZETAL(irun)        _RL VPSIH(irun),VZETAL(irun)
1592        real VZ0L(irun),VPSIH2(irun)        _RL VZ0L(irun),VPSIH2(irun)
1593        real VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)        _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
1594        real VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)        _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
1595        real VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)        _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
1596        real VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)        _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
1597        real VXNUM(irun),VDZETA1(irun),VDZETA2(irun)        _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
1598        real VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)        _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
1599        real VDPSIMC(irun),VDPSIHC(irun)        _RL VDPSIMC(irun),VDPSIHC(irun)
1600    
1601        real DZITRP(irun,nlev-1),STBFCN(irun,nlev)        _RL DZITRP(irun,nlev-1),STBFCN(irun,nlev)
1602        real XL0(irun,nlev),Q1(irun,nlev-1)        _RL XL0(irun,nlev),Q1(irun,nlev-1)
1603        real WRKIT1(irun,nlev-1)        _RL WRKIT1(irun,nlev-1)
1604        real WRKIT2(irun,nlev-1)        _RL WRKIT2(irun,nlev-1)
1605        real WRKIT3(irun,nlev-1)        _RL WRKIT3(irun,nlev-1)
1606        real WRKIT4(irun,nlev-1)        _RL WRKIT4(irun,nlev-1)
1607        INTEGER INT1(irun,nlev), INT2(irun,nlev-1)        INTEGER INT1(irun,nlev), INT2(irun,nlev-1)
1608    
1609        real vrt1con,pi,rsq2pi,p5sr,clh,vk,rvk,aitr,gbycp,fac1,fac2        _RL vrt1con,pi,rsq2pi,p5sr,clh,vk,rvk,aitr,gbycp,fac1,fac2
1610        real getcon,dum,errf        _RL getcon,dum,errf
1611        integer istnlv,nlevm1,nlevm2,nlevml,nlevp1,istnm1,istnm2,istnp1        integer istnlv,nlevm1,nlevm2,nlevml,nlevp1,istnm1,istnm2,istnp1
1612        integer istnml,istnmq,istlmq,nlevmq        integer istnml,istnmq,istlmq,nlevmq
1613        integer i,iter,init,n,nt,LL,L,Lp,Lp1,lmin,lminq,lminq1,ibit        integer i,iter,init,n,LL,L,Lp,Lp1,lmin,lminq,lminq1,ibit
1614    
1615        vk = getcon('VON KARMAN')        vk = getcon('VON KARMAN')
1616        rvk = 1./vk        rvk = 1./vk
# Line 1629  C Some space for variables to be used in Line 1623  C Some space for variables to be used in
1623        ISTNM2 = irun * NLEVM2        ISTNM2 = irun * NLEVM2
1624        ISTNP1 = irun * NLEVP1        ISTNP1 = irun * NLEVP1
1625        GBYCP  = GRAV / CP        GBYCP  = GRAV / CP
1626    
1627        VRT1CON = 1. + VIRTCON        VRT1CON = 1. + VIRTCON
1628        PI   = 4. * ATAN(1.)        PI   = 4. * ATAN(1.)
1629        RSQ2PI = 1./ ((2.*PI)**0.5)        RSQ2PI = 1./ ((2.*PI)**0.5)
# Line 1645  C -------------------------- Line 1639  C --------------------------
1639        IF(QBEG) INIT = 1        IF(QBEG) INIT = 1
1640  C SET DIAGNOSTIC LOGICALS AND INITIALIZE DIAGNOSTIC ARRAYS  C SET DIAGNOSTIC LOGICALS AND INITIALIZE DIAGNOSTIC ARRAYS
1641  C --------------------------------------------------------  C --------------------------------------------------------
1642        do I =1,istnlv  c     do I =1,istnlv
1643        wu(i,1) = 0.  c     wu(i,1) = 0.
1644        enddo  c     enddo
1645        do I =1,istnlv  c     do I =1,istnlv
1646        wv(i,1) = 0.  c     wv(i,1) = 0.
1647        enddo  c     enddo
1648        do I =1,istnlv  c     do I =1,istnlv
1649        eu(i,1) = 0.  c     eu(i,1) = 0.
1650        enddo  c     enddo
1651        do I =1,istnlv  c     do I =1,istnlv
1652        et(i,1) = 0.  c     et(i,1) = 0.
1653        enddo  c     enddo
1654    c     if (tprof) then
1655    c     DO I =1,ISTNM1
1656    c     XLDIAG(I,1) = 0.
1657    c     enddo
1658    c     endif
1659    c     do I =1,irun
1660    c     wu(i,nlev) = 0.
1661    c     enddo
1662    c     do I =1,irun
1663    c     wv(i,nlev) = 0.
1664    c     enddo
1665          DO L =1,NLEV
1666            DO I =1,irun
1667             wu(I,L) = 0.
1668             wv(I,L) = 0.
1669             eu(I,L) = 0.
1670             et(I,L) = 0.
1671            ENDDO
1672          ENDDO
1673        if (tprof) then        if (tprof) then
1674        DO I =1,ISTNM1         DO L =1,NLEVM1
1675        XLDIAG(I,1) = 0.          DO I =1,irun
1676        enddo           XLDIAG(I,L) = 0.
1677            ENDDO
1678           ENDDO
1679        endif        endif
1680        do I =1,irun        DO L =1,NLEVM1
1681        wu(i,nlev) = 0.          DO I =1,irun
1682        enddo           FREQDG(I,L) = 0.
1683        do I =1,irun          ENDDO
1684        wv(i,nlev) = 0.        ENDDO
       enddo  
1685        do I =1,irun        do I =1,irun
1686        scu(i) = 0.        scu(i) = 0.
1687        enddo        enddo
# Line 1683  C -------------------------------------- Line 1697  C --------------------------------------
1697        do I =1,irun        do I =1,irun
1698        sz0(i) = 0.        sz0(i) = 0.
1699        enddo        enddo
1700        do I =1,ISTNM1  c     do I =1,ISTNM1
1701        FREQDG(I,1) = 0.  c     FREQDG(I,1) = 0.
1702        enddo  c     enddo
1703        do I =1,irun        do I =1,irun
1704        stu2m(i) = 0.        stu2m(i) = 0.
1705        enddo        enddo
# Line 1712  C -------------------------------------- Line 1726  C --------------------------------------
1726        enddo        enddo
1727    
1728        IF (INIT.EQ.1) THEN        IF (INIT.EQ.1) THEN
1729         DO I = 1,ISTNM1  c      DO I = 1,ISTNM1
1730          XLSAVE(I,1) = 0.  c       XLSAVE(I,1) = 0.
1731          KHSAVE(I,1) = 0.  c       KHSAVE(I,1) = 0.
1732    c      ENDDO
1733           DO L =1,NLEVM1
1734            DO I =1,irun
1735            XLSAVE(I,L) = 0.
1736            KHSAVE(I,L) = 0.
1737            ENDDO
1738         ENDDO         ENDDO
1739         DO I = 1,irun         DO I = 1,irun
1740          CTSAVE(I) = 0.          CTSAVE(I) = 0.
# Line 1726  C -------------------------------------- Line 1746  C --------------------------------------
1746    
1747  C COMPUTE VERTICAL GRID  C COMPUTE VERTICAL GRID
1748  C ---------------------  C ---------------------
1749        DO 9038 I =1,ISTNLV  c     DO 9038 I =1,ISTNLV
1750        ADZ1(I,1) = (CP/GRAV)*(PLKE(I,2)-PLKE(I,1))  c     ADZ1(I,1) = (CP/GRAV)*(PLKE(I,2)-PLKE(I,1))
1751        ADZ1(I,1) = THV(I,1) * ADZ1(I,1)  c     ADZ1(I,1) = THV(I,1) * ADZ1(I,1)
1752        DZ1TMP(I,1) = ADZ1(I,1)  c     DZ1TMP(I,1) = ADZ1(I,1)
1753  9038  CONTINUE  9038  CONTINUE
1754        DO 9040 I =1,ISTNM1        DO L =1,NLEV
1755        ADZ2(I,1) = 0.5 * (ADZ1(I,1)+ADZ1(I,2))         DO I =1,irun
1756            ADZ1(I,L) = (CP/GRAV)*(PLKE(I,L+1)-PLKE(I,L))
1757            ADZ1(I,L) = THV(I,L) * ADZ1(I,L)
1758            DZ1TMP(I,L) = ADZ1(I,L)
1759           ENDDO
1760          ENDDO
1761    c     DO 9040 I =1,ISTNM1
1762    c     ADZ2(I,1) = 0.5 * (ADZ1(I,1)+ADZ1(I,2))
1763  9040  CONTINUE  9040  CONTINUE
1764          DO L =1,NLEVM1
1765           DO I =1,irun
1766            ADZ2(I,L) = 0.5 * (ADZ1(I,L)+ADZ1(I,L+1))
1767           ENDDO
1768          ENDDO
1769  C DEPTH HS OF SURFACE LAYER  C DEPTH HS OF SURFACE LAYER
1770  C -------------------------  C -------------------------
1771        DO 9042 I =1,irun        DO 9042 I =1,irun
# Line 1744  C ------------------------------- Line 1776  C -------------------------------
1776        DO 9044 I =1,irun        DO 9044 I =1,irun
1777        DZ3(I,1) = HALPHA * ADZ1(I,1)        DZ3(I,1) = HALPHA * ADZ1(I,1)
1778  9044  CONTINUE  9044  CONTINUE
1779        DO 9046 I =1,ISTNM2  c     DO 9046 I =1,ISTNM2
1780        DZ3(I,2) = ALPHA * ADZ1(I,2)  c     DZ3(I,2) = ALPHA * ADZ1(I,2)
1781  9046  CONTINUE  9046  CONTINUE
1782          DO L =2,NLEVM1
1783           DO I =1,irun
1784            DZ3(I,L) = ALPHA * ADZ1(I,L)
1785           ENDDO
1786          ENDDO
1787        DO 9048 I =1,irun        DO 9048 I =1,irun
1788        DZ3(I,NLEV) = ALPHA * HS(I)        DZ3(I,NLEV) = ALPHA * HS(I)
1789  9048  CONTINUE  9048  CONTINUE
1790    
1791  C VK * HEIGHTS AT MID AND EDGE LEVELS  C VK * HEIGHTS AT MID AND EDGE LEVELS
1792  C -----------------------------------  C -----------------------------------
1793        DO 9050 I =1,ISTNM1  c     DO 9050 I =1,ISTNM1
1794        TEMP(I,2) = VK * ADZ1(I,2)  c     TEMP(I,2) = VK * ADZ1(I,2)
1795  9050  CONTINUE  9050  CONTINUE
1796          DO L =2,NLEV
1797           DO I =1,irun
1798            TEMP(I,L) = VK * ADZ1(I,L)
1799           ENDDO
1800          ENDDO
1801        DO 9052 I =1,irun        DO 9052 I =1,irun
1802        VKZE(I,NLEVM1) = TEMP(I,NLEV)        VKZE(I,NLEVM1) = TEMP(I,NLEV)
1803  9052  CONTINUE  9052  CONTINUE
# Line 1766  C ----------------------------------- Line 1808  C -----------------------------------
1808        VKZE(I,L) = VKZE(I,LP1) + TEMP(I,LP1)        VKZE(I,L) = VKZE(I,LP1) + TEMP(I,LP1)
1809  9054  CONTINUE  9054  CONTINUE
1810   100  CONTINUE   100  CONTINUE
1811        DO 9056 I =1,ISTNM1  c     DO 9056 I =1,ISTNM1
1812        VKZM(I,1) = VKZE(I,1) - 0.5 * TEMP(I,2)  c     VKZM(I,1) = VKZE(I,1) - 0.5 * TEMP(I,2)
1813  9056  CONTINUE  9056  CONTINUE
1814          DO L =1,NLEVM1
1815           DO I =1,irun
1816            VKZM(I,L) = VKZE(I,L) - 0.5 * TEMP(I,L+1)
1817           ENDDO
1818          ENDDO
1819  C COMPUTE RHO BY DZ AT MID AND EDGE LEVELS  C COMPUTE RHO BY DZ AT MID AND EDGE LEVELS
1820  C ----------------------------------------  C ----------------------------------------
1821        DO 200 L = 1,NLEVM1        DO 200 L = 1,NLEVM1
# Line 1780  C -------------------------------------- Line 1827  C --------------------------------------
1827         RHODZ2(I,L) = RHODZ2(I,L) + FAC2 * THV(I,L)         RHODZ2(I,L) = RHODZ2(I,L) + FAC2 * THV(I,L)
1828  9058  CONTINUE  9058  CONTINUE
1829   200  CONTINUE   200  CONTINUE
1830        DO 9060 I =1,ISTNM1  c     DO 9060 I =1,ISTNM1
1831         RHODZ2(I,1) = (RGAS*0.01) * RHODZ2(I,1)  c      RHODZ2(I,1) = (RGAS*0.01) * RHODZ2(I,1)
1832         TEMP(I,1)   = PLKE(I,2) * ADZ2(I,1)  c      TEMP(I,1)   = PLKE(I,2) * ADZ2(I,1)
1833         RHODZ2(I,1) = TEMP(I,1) * RHODZ2(I,1)  c      RHODZ2(I,1) = TEMP(I,1) * RHODZ2(I,1)
1834         RHODZ2(I,1) = PLE(I,2) / RHODZ2(I,1)  c      RHODZ2(I,1) = PLE(I,2) / RHODZ2(I,1)
1835         RHOZPK(I,1) = RHODZ2(I,1) * PLKE(I,2)  c      RHOZPK(I,1) = RHODZ2(I,1) * PLKE(I,2)
1836         RHODZ1(I,1) = (RGAS*0.01) * THV(I,2)  c      RHODZ1(I,1) = (RGAS*0.01) * THV(I,2)
1837         TEMP(I,1)   = PLK(I,2) * ADZ1(I,2)  c      TEMP(I,1)   = PLK(I,2) * ADZ1(I,2)
1838         RHODZ1(I,1) = TEMP(I,1) * RHODZ1(I,1)  c      RHODZ1(I,1) = TEMP(I,1) * RHODZ1(I,1)
1839         RHODZ1(I,1) = PL(I,2) / RHODZ1(I,1)  c      RHODZ1(I,1) = PL(I,2) / RHODZ1(I,1)
1840  9060  CONTINUE  9060  CONTINUE
1841          DO L =1,NLEVM1
1842           DO I =1,irun
1843            RHODZ2(I,L) = (RGAS*0.01) * RHODZ2(I,L)
1844            TEMP(I,L)   = PLKE(I,L+1) * ADZ2(I,L)
1845            RHODZ2(I,L) = TEMP(I,L) * RHODZ2(I,L)
1846            RHODZ2(I,L) = PLE(I,L+1) / RHODZ2(I,L)
1847            RHOZPK(I,L) = RHODZ2(I,L) * PLKE(I,L+1)
1848            RHODZ1(I,L) = (RGAS*0.01) * THV(I,L+1)
1849            TEMP(I,L)   = PLK(I,L+1) * ADZ1(I,L+1)
1850            RHODZ1(I,L) = TEMP(I,L) * RHODZ1(I,L)
1851            RHODZ1(I,L) = PL(I,L+1) / RHODZ1(I,L)
1852           ENDDO
1853          ENDDO
1854  C COMPUTE FLXFAC FOR LAYERS AND EDGES  C COMPUTE FLXFAC FOR LAYERS AND EDGES
1855  C COMPUTE DTG / DT DUE TO RADIATION AND HEAT CONDUCTION THROUGH ICE  C COMPUTE DTG / DT DUE TO RADIATION AND HEAT CONDUCTION THROUGH ICE
1856  C -----------------------------------------------------------------  C -----------------------------------------------------------------
1857        DO 9062 I =1,ISTNLV  c     DO 9062 I =1,ISTNLV
1858        FLXFPK(I,1) = PLE(I,2) - PLE(I,1)  c     FLXFPK(I,1) = PLE(I,2) - PLE(I,1)
1859        FLXFPK(I,1) = FLXFPK(I,1) * PLK(I,1)  c     FLXFPK(I,1) = FLXFPK(I,1) * PLK(I,1)
1860        FLXFPK(I,1) = (GRAV*DTAU*0.01) / FLXFPK(I,1)  c     FLXFPK(I,1) = (GRAV*DTAU*0.01) / FLXFPK(I,1)
1861  9062  CONTINUE  9062  CONTINUE
1862          DO L =1,NLEV
1863           DO I =1,irun
1864            FLXFPK(I,L) = PLE(I,L+1) - PLE(I,L)
1865            FLXFPK(I,L) = FLXFPK(I,L) * PLK(I,L)
1866            FLXFPK(I,L) = (GRAV*DTAU*0.01) / FLXFPK(I,L)
1867           ENDDO
1868          ENDDO
1869        DO 9064 I =1,irun        DO 9064 I =1,irun
1870        FLXFPK(I,NLEVP1) = 0.        FLXFPK(I,NLEVP1) = 0.
1871  9064  CONTINUE  9064  CONTINUE
1872        DO 9066 I =1,irun        DO 9066 I =1,irun
1873        IF (IWATER(I).EQ.0 ) FLXFPK(I,NLEVP1) = 1. / PLKE(I,NLEVP1)        IF (IWATER(I).EQ.0 ) FLXFPK(I,NLEVP1) = 1. / PLKE(I,NLEVP1)
1874  9066  CONTINUE  9066  CONTINUE
1875        DO 9068 I =1,ISTNLV  c     DO 9068 I =1,ISTNLV
1876        FLXFAC(I,1) = FLXFPK(I,1) * PLK(I,1)  c     FLXFAC(I,1) = FLXFPK(I,1) * PLK(I,1)
1877  9068  CONTINUE  9068  CONTINUE
1878          DO L =1,NLEV
1879           DO I =1,irun
1880            FLXFAC(I,L) = FLXFPK(I,L) * PLK(I,L)
1881           ENDDO
1882          ENDDO
1883        DO 9070 I =1,irun        DO 9070 I =1,irun
1884        FLXFAC(I,NLEVP1) = FLXFPK(I,NLEVP1)        FLXFAC(I,NLEVP1) = FLXFPK(I,NLEVP1)
1885  9070  CONTINUE  9070  CONTINUE
1886        DO 9074 I =1,irun        DO 9074 I =1,irun
1887        FLXFPK(I,NLEVP1) = CP * FLXFPK(I,NLEVP1)        FLXFPK(I,NLEVP1) = CP * FLXFPK(I,NLEVP1)
1888  9074  CONTINUE  9074  CONTINUE
1889        DO 9076 I =1,ISTNM1  c     DO 9076 I =1,ISTNM1
1890        FLXFCE(I,1) = PL(I,2) - PL(I,1)  c     FLXFCE(I,1) = PL(I,2) - PL(I,1)
1891  9076  CONTINUE  9076  CONTINUE
1892        DO 9078 I =1,ISTNM1  c     DO 9078 I =1,ISTNM1
1893        FLXFCE(I,1) = (GRAV*DTAU*0.01) / FLXFCE(I,1)  c     FLXFCE(I,1) = (GRAV*DTAU*0.01) / FLXFCE(I,1)
1894  9078  CONTINUE  9078  CONTINUE
1895          DO L =1,NLEVM1
1896           DO I =1,irun
1897            FLXFCE(I,L) = PL(I,L+1) - PL(I,L)
1898            FLXFCE(I,L) = (GRAV*DTAU*0.01) / FLXFCE(I,L)
1899           ENDDO
1900          ENDDO
1901  C COMPUTE RECIPROCALS OF DZ1, DZ2, HS  C COMPUTE RECIPROCALS OF DZ1, DZ2, HS
1902  C -----------------------------------  C -----------------------------------
1903        DO 9084 I =1,ISTNLV  c     DO 9084 I =1,ISTNLV
1904        ADZ1(I,1) = 1. / ADZ1(I,1)  c     ADZ1(I,1) = 1. / ADZ1(I,1)
1905  9084  CONTINUE  9084  CONTINUE
1906        DO 9086 I =1,ISTNM1        DO L =1,NLEV
1907        ADZ2(I,1) = 1. / ADZ2(I,1)         DO I =1,irun
1908            ADZ1(I,L) = 1. / ADZ1(I,L)
1909           ENDDO
1910          ENDDO
1911    c     DO 9086 I =1,ISTNM1
1912    c     ADZ2(I,1) = 1. / ADZ2(I,1)
1913  9086  CONTINUE  9086  CONTINUE
1914          DO L =1,NLEVM1
1915           DO I =1,irun
1916            ADZ2(I,L) = 1. / ADZ2(I,L)
1917           ENDDO
1918          ENDDO
1919        DO 9088 I =1,irun        DO 9088 I =1,irun
1920        AHS(I) = 1. / HS(I)        AHS(I) = 1. / HS(I)
1921  9088  CONTINUE  9088  CONTINUE
1922  C COMPUTE GRADIENTS OF P**KAPPA  C COMPUTE GRADIENTS OF P**KAPPA
1923  C -----------------------------  C -----------------------------
1924        DO 9090 I =1,ISTNM1  c     DO 9090 I =1,ISTNM1
1925        DPK(I,1) = ( PLK(I,2)-PLK(I,1) ) * ADZ2(I,1)  c     DPK(I,1) = ( PLK(I,2)-PLK(I,1) ) * ADZ2(I,1)
1926  9090  CONTINUE  9090  CONTINUE
1927          DO L =1,NLEVM1
1928           DO I =1,irun
1929            DPK(I,L) = ( PLK(I,L+1)-PLK(I,L) ) * ADZ2(I,L)
1930           ENDDO
1931          ENDDO
1932        DO 9092 I =1,irun        DO 9092 I =1,irun
1933        DPK(I,NLEV) = GBYCP / THV(I,NLEV)        DPK(I,NLEV) = GBYCP / THV(I,NLEV)
1934  9092  CONTINUE  9092  CONTINUE
1935  C INITIALIZE Q ARRAY  C INITIALIZE Q ARRAY
1936  C ------------------  C ------------------
1937        DO 9094 I =1,ISTNM1  c     DO 9094 I =1,ISTNM1
1938        Q(I,1) = 2. * QQ(I,1)  c     Q(I,1) = 2. * QQ(I,1)
1939        Q(I,1) = SQRT( Q(I,1) )  c     Q(I,1) = SQRT( Q(I,1) )
1940  9094  CONTINUE  9094  CONTINUE
1941          DO L =1,NLEVM1
1942           DO I =1,irun
1943            Q(I,L) = 2. * QQ(I,L)
1944            Q(I,L) = SQRT( Q(I,L) )
1945           ENDDO
1946          ENDDO
1947        FIRST = .TRUE.        FIRST = .TRUE.
1948        LAST = .FALSE.        LAST = .FALSE.
1949  C**********************************************************************  C**********************************************************************
# Line 1943  CE Line 2042  CE
2042  C  C
2043  C GRADIENTS  C GRADIENTS
2044  C ---------  C ---------
2045        DO 9100 I =1,ISTNM1  c     DO 9100 I =1,ISTNM1
2046          DU(I,1) = (   U(I,1)-  U(I,2) ) * ADZ2(I,1)  c       DU(I,1) = (   U(I,1)-  U(I,2) ) * ADZ2(I,1)
2047          DV(I,1) = (   V(I,1)-  V(I,2) ) * ADZ2(I,1)  c       DV(I,1) = (   V(I,1)-  V(I,2) ) * ADZ2(I,1)
2048  9100  CONTINUE  9100  CONTINUE
2049          DO L =1,NLEVM1
2050           DO I =1,irun
2051            DU(I,L) = (   U(I,L)-  U(I,L+1) ) * ADZ2(I,L)
2052            DV(I,L) = (   V(I,L)-  V(I,L+1) ) * ADZ2(I,L)
2053           ENDDO
2054          ENDDO
2055    
2056    
2057  C    NEW CODE FOR MOIST BOUNDARY LAYER - NEW CALCULATION OF DTHV  C    NEW CODE FOR MOIST BOUNDARY LAYER - NEW CALCULATION OF DTHV
2058  C  C
2059        IF(ITER.EQ.1) THEN        IF(ITER.EQ.1) THEN
2060         DO I = 1,ISTNM1  c      DO I = 1,ISTNM1
2061          XL(I,1) = XLSAVE(I,1)  c       XL(I,1) = XLSAVE(I,1)
2062    c      ENDDO
2063           DO L =1,NLEVM1
2064            DO I =1,irun
2065             XL(I,L) = XLSAVE(I,L)
2066            ENDDO
2067         ENDDO         ENDDO
2068        ENDIF        ENDIF
2069  C  C
2070        DO I =1,ISTNM1  c     DO I =1,ISTNM1
2071         DTH(I,1) = ( TH(I,1)-TH(I,2) ) * ADZ2(I,1)  c      DTH(I,1) = ( TH(I,1)-TH(I,2) ) * ADZ2(I,1)
2072         DSH(I,1) = ( SH(I,1)-SH(I,2) ) * ADZ2(I,1)  c      DSH(I,1) = ( SH(I,1)-SH(I,2) ) * ADZ2(I,1)
2073         TL(I,1) = TH(I,1)*PLK(I,1)  c      TL(I,1) = TH(I,1)*PLK(I,1)
2074    c     ENDDO
2075          DO L =1,NLEVM1
2076           DO I =1,irun
2077            DTH(I,L) = ( TH(I,L)-TH(I,L+1) ) * ADZ2(I,L)
2078            DSH(I,L) = ( SH(I,L)-SH(I,L+1) ) * ADZ2(I,L)
2079            TL(I,L) = TH(I,L)*PLK(I,L)
2080           ENDDO
2081        ENDDO        ENDDO
2082        DO LL = 1,NLEVM1          DO LL = 1,NLEVM1
2083        DO I = 1,irun          DO I = 1,irun
2084         call qsat ( tl(i,LL),pl(i,LL),shsat(i,LL),dum,.false. )         call qsat ( tl(i,LL),pl(i,LL),shsat(i,LL),dum,.false. )
2085        ENDDO        ENDDO
2086        ENDDO        ENDDO
2087        DO I = 1,ISTNM1  c     DO I = 1,ISTNM1
2088         BB(I,1) = FACEPS*SHSAT(I,1)/(TL(I,1)*TL(I,1))  c      BB(I,1) = FACEPS*SHSAT(I,1)/(TL(I,1)*TL(I,1))
2089         AA(I,1) = 1. / (1. + CLH * BB(I,1) )  c      AA(I,1) = 1. / (1. + CLH * BB(I,1) )
2090  COMMM  BB(I,1) = BB(I,1) * AA(I,1) * plke(I,2)  cCOMMM BB(I,1) = BB(I,1) * AA(I,1) * plke(I,2)
2091         BB(I,1) = BB(I,1) * AA(I,1)  c      BB(I,1) = BB(I,1) * AA(I,1)
2092         SBAR(I,1) = AA(I,1) * (SH(I,1) - SHSAT(I,1))  c      SBAR(I,1) = AA(I,1) * (SH(I,1) - SHSAT(I,1))
2093    c     ENDDO
2094          DO L =1,NLEVM1
2095           DO I =1,irun
2096            BB(I,L) = FACEPS*SHSAT(I,L)/(TL(I,L)*TL(I,L))
2097            AA(I,L) = 1. / (1. + CLH * BB(I,L) )
2098    COMMM   BB(I,L) = BB(I,L) * AA(I,L) * plke(I,L+1)
2099            BB(I,L) = BB(I,L) * AA(I,L)
2100            SBAR(I,L) = AA(I,L) * (SH(I,L) - SHSAT(I,L))
2101           ENDDO
2102        ENDDO        ENDDO
2103        DO I = 1,irun        DO I = 1,irun
2104  COMMM  SSDEV(I,1) = XL(I,1)*(AA(I,1)*DSH(I,1)-BB(I,1)*DTH(I,1))  COMMM  SSDEV(I,1) = XL(I,1)*(AA(I,1)*DSH(I,1)-BB(I,1)*DTH(I,1))
# Line 1982  COMMM  SSDEV(I,1) = XL(I,1)*(AA(I,1)*DSH Line 2108  COMMM  SSDEV(I,1) = XL(I,1)*(AA(I,1)*DSH
2108         SVAR(I,1) = SQRT(SSDEV(I,1))         SVAR(I,1) = SQRT(SSDEV(I,1))
2109         IF ( SVAR(I,1).LT.Z1PEM25) SVAR(I,1) = Z1PEM25         IF ( SVAR(I,1).LT.Z1PEM25) SVAR(I,1) = Z1PEM25
2110        ENDDO        ENDDO
2111        DO I = 1,ISTNM2  c     DO I = 1,ISTNM2
2112  COMMM  SSDEV(I,2) = XL(I,1)*(AA(I,2)*DSH(I,1)-BB(I,2)*DTH(I,1))  cCOMMM SSDEV(I,2) = XL(I,1)*(AA(I,2)*DSH(I,1)-BB(I,2)*DTH(I,1))
2113         SSDEV(I,2) = XL(I,1)*(AA(I,2)*DSH(I,1)-  c      SSDEV(I,2) = XL(I,1)*(AA(I,2)*DSH(I,1)-
2114       1              BB(I,2)*plke(I,2)*DTH(I,1))  c    1              BB(I,2)*plke(I,2)*DTH(I,1))
2115         SSDEV(I,2) = B2 * KHSAVE(I,1) * SSDEV(I,1) * SSDEV(I,1)  c      SSDEV(I,2) = B2 * KHSAVE(I,1) * SSDEV(I,1) * SSDEV(I,1)
2116         SVAR(I,2) = SQRT(SSDEV(I,2))  c      SVAR(I,2) = SQRT(SSDEV(I,2))
2117  COMMM  SSDEV(I,2) = XL(I,2)*(AA(I,2)*DSH(I,2)-BB(I,2)*DTH(I,2))  cCOMMM SSDEV(I,2) = XL(I,2)*(AA(I,2)*DSH(I,2)-BB(I,2)*DTH(I,2))
2118         SSDEV(I,2) = XL(I,2)*(AA(I,2)*DSH(I,2)-  c      SSDEV(I,2) = XL(I,2)*(AA(I,2)*DSH(I,2)-
2119       1              BB(I,2)*plke(I,3)*DTH(I,2))  c    1              BB(I,2)*plke(I,3)*DTH(I,2))
2120         SSDEV(I,2) = B2 * KHSAVE(I,2) * SSDEV(I,2) * SSDEV(I,2)  c      SSDEV(I,2) = B2 * KHSAVE(I,2) * SSDEV(I,2) * SSDEV(I,2)
2121         TEMP(I,2) = SQRT(SSDEV(I,2))  c      TEMP(I,2) = SQRT(SSDEV(I,2))
2122         SVAR(I,2) = (1./2.) * (SVAR(I,2) + TEMP(I,2))  c      SVAR(I,2) = (1./2.) * (SVAR(I,2) + TEMP(I,2))
2123         IF ( SVAR(I,2).LT.Z1PEM25) SVAR(I,2) = Z1PEM25  c      IF ( SVAR(I,2).LT.Z1PEM25) SVAR(I,2) = Z1PEM25
2124        ENDDO  c     ENDDO
2125        DO I = 1,ISTNM1        DO L =2,NLEVM1
2126         Q1M(I,1) = SBAR(I,1) / SVAR(I,1)         DO I =1,irun
2127         FCC(I,1) = (1./2.) * ( 1. + ERRF( P5SR*Q1M(I,1) ) )  COMMM   SSDEV(I,L) = XL(I,L-1)*(AA(I,L)*DSH(I,L-1)-BB(I,L)*DTH(I,L-1))
2128         SHL(I,1) =  FCC(I,1) * SBAR(I,1)          SSDEV(I,L) = XL(I,L-1)*(AA(I,L)*DSH(I,L-1)-
2129         ARG(I,1)  = (1./2.)*Q1M(I,1)*Q1M(I,1)       1               BB(I,L)*plke(I,L)*DTH(I,L-1))
2130         IF(ARG(I,1).LE.ARGMAX)          SSDEV(I,L) = B2 * KHSAVE(I,L-1) * SSDEV(I,L-1) * SSDEV(I,L-1)
2131       1       SHL(I,1) =  SHL(I,1)+RSQ2PI*SVAR(I,1)*EXP(-ARG(I,1))          SVAR(I,L) = SQRT(SSDEV(I,L))
2132         BETAT(I,1) = 1. + VIRTCON * SH(I,1) - VRT1CON * SHL(I,1)  COMMM   SSDEV(I,L) = XL(I,L)*(AA(I,L)*DSH(I,L)-BB(I,L)*DTH(I,L))
2133         BETAW(I,1) = VIRTCON *          SSDEV(I,L) = XL(I,L)*(AA(I,L)*DSH(I,L)-
2134       1               ( TH(I,1) + (CLH/plk(I,1)) * SHL(I,1) )       1               BB(I,L)*plke(I,L+1)*DTH(I,L))
2135         BETAL(I,1) = ( 1. + VIRTCON*SH(I,1) - TWO*VRT1CON*SHL(I,1) )          SSDEV(I,L) = B2 * KHSAVE(I,L) * SSDEV(I,L) * SSDEV(I,L)
2136       1    * (CLH/plke(I,2)) - VRT1CON * TH(I,1)          TEMP(I,L) = SQRT(SSDEV(I,L))
2137  COMMM  BETAT1(I,1) = BETAT(I,1) - BB(I,1) * FCC(I,1) * BETAL(I,1)          SVAR(I,L) = (1./2.) * (SVAR(I,L) + TEMP(I,L))
2138         BETAT1(I,1) = BETAT(I,1) -          IF ( SVAR(I,L).LT.Z1PEM25) SVAR(I,L) = Z1PEM25
2139       1              BB(I,1)*plk(i,1) * FCC(I,1) * BETAL(I,1)         ENDDO
2140         BETAW1(I,1) = BETAW(I,1) + AA(I,1) * FCC(I,1) * BETAL(I,1)        ENDDO
2141        ENDDO  c     DO I = 1,ISTNM1
2142        DO I = 1,ISTNM1  c      Q1M(I,1) = SBAR(I,1) / SVAR(I,1)
2143         DTHV(I,1) =  (1./2.)*((BETAT1(I,1)+BETAT1(I,2))*DTH(I,1)  c      FCC(I,1) = (1./2.) * ( 1. + ERRF( P5SR*Q1M(I,1) ) )
2144       1              + (BETAW1(I,1)+BETAW1(I,2))*DSH(I,1))  c      SHL(I,1) =  FCC(I,1) * SBAR(I,1)
2145    c      ARG(I,1)  = (1./2.)*Q1M(I,1)*Q1M(I,1)
2146    c      IF(ARG(I,1).LE.ARGMAX)
2147    c    1       SHL(I,1) =  SHL(I,1)+RSQ2PI*SVAR(I,1)*EXP(-ARG(I,1))
2148    c      BETAT(I,1) = 1. + VIRTCON * SH(I,1) - VRT1CON * SHL(I,1)
2149    c      BETAW(I,1) = VIRTCON *
2150    c    1               ( TH(I,1) + (CLH/plk(I,1)) * SHL(I,1) )
2151    c      BETAL(I,1) = ( 1. + VIRTCON*SH(I,1) - TWO*VRT1CON*SHL(I,1) )
2152    c    1    * (CLH/plke(I,2)) - VRT1CON * TH(I,1)
2153    cCOMMM BETAT1(I,1) = BETAT(I,1) - BB(I,1) * FCC(I,1) * BETAL(I,1)
2154    c      BETAT1(I,1) = BETAT(I,1) -
2155    c    1              BB(I,1)*plk(i,1) * FCC(I,1) * BETAL(I,1)
2156    c      BETAW1(I,1) = BETAW(I,1) + AA(I,1) * FCC(I,1) * BETAL(I,1)
2157    c     ENDDO
2158          DO L =1,NLEVM1
2159           DO I =1,irun
2160            Q1M(I,L) = SBAR(I,L) / SVAR(I,L)
2161            FCC(I,L) = (1./2.) * ( 1. + ERRF( P5SR*Q1M(I,L) ) )
2162            SHL(I,L) =  FCC(I,L) * SBAR(I,L)
2163            ARG(I,L)  = (1./2.)*Q1M(I,L)*Q1M(I,L)
2164            IF(ARG(I,L).LE.ARGMAX)
2165         1        SHL(I,L) =  SHL(I,L)+RSQ2PI*SVAR(I,L)*EXP(-ARG(I,L))
2166            BETAT(I,L) = 1. + VIRTCON * SH(I,L) - VRT1CON * SHL(I,L)
2167            BETAW(I,L) = VIRTCON *
2168         1                ( TH(I,L) + (CLH/plk(I,L)) * SHL(I,L) )
2169            BETAL(I,L) = ( 1. + VIRTCON*SH(I,L) - TWO*VRT1CON*SHL(I,L) )
2170         1     * (CLH/plke(I,L+1)) - VRT1CON * TH(I,L)
2171    COMMM   BETAT1(I,L) = BETAT(I,L) - BB(I,L) * FCC(I,L) * BETAL(I,L)
2172            BETAT1(I,L) = BETAT(I,L) -
2173         1               BB(I,L)*plk(i,L) * FCC(I,L) * BETAL(I,L)
2174            BETAW1(I,L) = BETAW(I,L) + AA(I,L) * FCC(I,L) * BETAL(I,L)
2175           ENDDO
2176          ENDDO
2177    c     DO I = 1,ISTNM1
2178    c      DTHV(I,1) =  (1./2.)*((BETAT1(I,1)+BETAT1(I,2))*DTH(I,1)
2179    c    1              + (BETAW1(I,1)+BETAW1(I,2))*DSH(I,1))
2180    c     ENDDO
2181          DO L =1,NLEVM1
2182           DO I =1,irun
2183            DTHV(I,L) =  (1./2.)*((BETAT1(I,L)+BETAT1(I,L+1))*DTH(I,L)
2184         1               + (BETAW1(I,L)+BETAW1(I,L+1))*DSH(I,L))
2185           ENDDO
2186        ENDDO        ENDDO
2187    
2188  C GRADIENTS AT THE TOP OF THE SURFACE LAYER  C GRADIENTS AT THE TOP OF THE SURFACE LAYER
2189  C -----------------------------------------  C -----------------------------------------
2190        DO 9102 I =1,irun        DO 9102 I =1,irun
# Line 2030  C -------------------------------------- Line 2197  C --------------------------------------
2197    
2198  C CALCULATE BRUNT-VAISALA FREQUENCIES, SHEARS, RICHARDSON NUMBERS  C CALCULATE BRUNT-VAISALA FREQUENCIES, SHEARS, RICHARDSON NUMBERS
2199  C ---------------------------------------------------------------  C ---------------------------------------------------------------
2200        DO 9104 I =1,ISTNLV  c     DO 9104 I =1,ISTNLV
2201        STRT(I,1) = CP * DTHV(I,1) * DPK(I,1)  c     STRT(I,1) = CP * DTHV(I,1) * DPK(I,1)
2202        DW2(I,1) = DU(I,1) * DU(I,1) + DV(I,1) * DV(I,1)  c     DW2(I,1) = DU(I,1) * DU(I,1) + DV(I,1) * DV(I,1)
2203        IF ( DW2(I,1) .LE. 1.e-4 ) DW2(I,1) = 1.e-4  c     IF ( DW2(I,1) .LE. 1.e-4 ) DW2(I,1) = 1.e-4
2204        RI(I,1) = STRT(I,1) / DW2(I,1)  c     RI(I,1) = STRT(I,1) / DW2(I,1)
2205  9104  CONTINUE  9104  CONTINUE
2206          DO L =1,NLEV
2207           DO I =1,irun
2208            STRT(I,L) = CP * DTHV(I,L) * DPK(I,L)
2209            DW2(I,L) = DU(I,L) * DU(I,L) + DV(I,L) * DV(I,L)
2210            IF ( DW2(I,L) .LE. 1.e-4 ) DW2(I,L) = 1.e-4
2211            RI(I,L) = STRT(I,L) / DW2(I,L)
2212           ENDDO
2213          ENDDO
2214  C FILL RICHARDSON NUMBER AND SURFACE WIND DIAGNOSTICS  C FILL RICHARDSON NUMBER AND SURFACE WIND DIAGNOSTICS
2215  C       (THOSE NEEDED FROM FIRST TRBFLX ITERATION)  C       (THOSE NEEDED FROM FIRST TRBFLX ITERATION)
2216  C ---------------------------------------------------  C ---------------------------------------------------
2217        DO 9106 I =1,ISTNM1  c     DO 9106 I =1,ISTNM1
2218        SRI(I,1) = RI(I,1)  c     SRI(I,1) = RI(I,1)
2219  9106  CONTINUE  9106  CONTINUE
2220          DO L =1,NLEVM1
2221           DO I =1,irun
2222            SRI(I,L) = RI(I,L)
2223           ENDDO
2224          ENDDO
2225        DO 9108 I =1,irun        DO 9108 I =1,irun
2226        SRI(I,NLEV) = RIB(I)        SRI(I,NLEV) = RIB(I)
2227  9108  CONTINUE  9108  CONTINUE
# Line 2050  C -------------------------------------- Line 2230  C --------------------------------------
2230  9110  CONTINUE  9110  CONTINUE
2231  C INITIALIZE KH, KM, QE AND P3 AND ELIMINATE SMALL QQ  C INITIALIZE KH, KM, QE AND P3 AND ELIMINATE SMALL QQ
2232  C ---------------------------------------------------  C ---------------------------------------------------
2233        DO 9112 I =1,ISTNM1  c     DO 9112 I =1,ISTNM1
2234        KH(I,1) = 0.  c     KH(I,1) = 0.
2235        KM(I,1) = 0.  c     KM(I,1) = 0.
2236        QQE(I,1) = 0.  c     QQE(I,1) = 0.
2237        QE(I,1) = 0.  c     QE(I,1) = 0.
2238        P3(I,1) = 0.  c     P3(I,1) = 0.
2239  9112  CONTINUE  9112  CONTINUE
2240        DO 9414 I = 1,ISTNM1  c     DO 9414 I = 1,ISTNM1
2241         IBITSTB(I,1) = 0  c      IBITSTB(I,1) = 0
2242   9414 CONTINUE  c9414 CONTINUE
2243        DO 9314 I = 1,ISTNM1  c     DO 9314 I = 1,ISTNM1
2244         IF ( QQ(I,1) .GT. 1.e-8 ) THEN  c      IF ( QQ(I,1) .GT. 1.e-8 ) THEN
2245          INTQ(I,1) = 1  c       INTQ(I,1) = 1
2246         ELSE  c      ELSE
2247          INTQ(I,1) = 0  c       INTQ(I,1) = 0
2248         ENDIF  c      ENDIF
2249   9314 CONTINUE  c9314 CONTINUE
2250        DO 9114 I = 1,ISTNM1  c     DO 9114 I = 1,ISTNM1
2251         IF ( QQ(I,1).LE.1.e-8 ) THEN  c      IF ( QQ(I,1).LE.1.e-8 ) THEN
2252          QQ(I,1) = 0.  c       QQ(I,1) = 0.
2253          Q(I,1) = 0.  c       Q(I,1) = 0.
2254         ENDIF  c      ENDIF
2255   9114 CONTINUE  c9114 CONTINUE
2256          DO L =1,NLEVM1
2257           DO I =1,irun
2258            KH(I,L) = 0.
2259            KM(I,L) = 0.
2260            QQE(I,L) = 0.
2261            QE(I,L) = 0.
2262            P3(I,L) = 0.
2263            IBITSTB(I,L) = 0
2264            IF ( QQ(I,L) .GT. 1.e-8 ) THEN
2265             INTQ(I,L) = 1
2266            ELSE
2267             INTQ(I,L) = 0
2268            ENDIF
2269            IF ( QQ(I,L).LE.1.e-8 ) THEN
2270             QQ(I,L) = 0.
2271             Q(I,L) = 0.
2272            ENDIF
2273           ENDDO
2274          ENDDO
2275  C  C
2276         DO 300 LMINQ = 1,NLEVM1         DO 300 LMINQ = 1,NLEVM1
2277          IBIT = 0          IBIT = 0
# Line 2102  C -------------------------------------- Line 2301  C --------------------------------------
2301  C FOR INITIAL START ONLY : USE EQUILIBRIUM MODEL  C FOR INITIAL START ONLY : USE EQUILIBRIUM MODEL
2302  C ----------------------------------------------  C ----------------------------------------------
2303        IF ( INIT .EQ. 1 ) THEN        IF ( INIT .EQ. 1 ) THEN
2304        DO 9180 I =1,ISTNM1  c     DO 9180 I =1,ISTNM1
2305         QQ(I,1) = QQE(I,1)  c      QQ(I,1) = QQE(I,1)
2306         Q(I,1) = QE(I,1)  c      Q(I,1) = QE(I,1)
2307  9180  CONTINUE  9180  CONTINUE
2308          DO L =1,NLEVM1
2309           DO I =1,irun
2310            QQ(I,L) = QQE(I,L)
2311            Q(I,L) = QE(I,L)
2312           ENDDO
2313          ENDDO
2314         INIT = 2         INIT = 2
2315          CALL TRBLEN(STRT,DW2,DZ3,Q,VKZE,VKZM,DTHV,DPK,DU,DV,XL,QXLM,          CALL TRBLEN(STRT,DW2,DZ3,Q,VKZE,VKZM,DTHV,DPK,DU,DV,XL,QXLM,
2316       1    NLEV,INIT,LMIN,LMINQ,LMINQ1,CP,INT1,INT2,       1    NLEV,INIT,LMIN,LMINQ,LMINQ1,CP,INT1,INT2,
# Line 2117  C DIMENSIONLESS COEFFS AND P3 (Q LE QE) Line 2322  C DIMENSIONLESS COEFFS AND P3 (Q LE QE)
2322  C -------------------------------------  C -------------------------------------
2323        IF( LMIN .LT. NLEV ) THEN        IF( LMIN .LT. NLEV ) THEN
2324         ISTNML = irun * NLEVML         ISTNML = irun * NLEVML
2325         DO 9320 I = 1,ISTNML  c      DO 9320 I = 1,ISTNML
2326          IF (  (IBITSTB(I,LMIN).EQ.1)  .AND.  c       IF (  (IBITSTB(I,LMIN).EQ.1)  .AND.
2327       1            ( Q(I,LMIN) .LE. QE(I,LMIN) ) ) THEN  c    1            ( Q(I,LMIN) .LE. QE(I,LMIN) ) ) THEN
2328           IBITSTB(I,LMIN) = 1  c        IBITSTB(I,LMIN) = 1
2329          ELSE  c       ELSE
2330           IBITSTB(I,LMIN) = 0  c        IBITSTB(I,LMIN) = 0
2331          ENDIF  c       ENDIF
2332   9320  CONTINUE  c9320  CONTINUE
2333         DO 9220 I = 1,ISTNML         DO L =LMIN,NLEVM1
2334          IF(IBITSTB(I,LMIN).EQ.1 ) THEN          DO I =1,irun
2335           TEMP(I,LMIN) = Q(I,LMIN) / QE(I,LMIN)           IF (  (IBITSTB(I,L).EQ.1)  .AND.
2336           KH(I,LMIN) = TEMP(I,LMIN) * KH(I,LMIN)       1             ( Q(I,L) .LE. QE(I,L) ) ) THEN
2337           KM(I,LMIN) = TEMP(I,LMIN) * KM(I,LMIN)            IBITSTB(I,L) = 1
2338          ENDIF           ELSE
2339          TEMP(I,LMIN) = 0.01 * QQE(I,LMIN)            IBITSTB(I,L) = 0
2340          IF((IBITSTB(I,LMIN).EQ.1) .AND.           ENDIF
2341       1        ( QQ(I,LMIN) .LE. TEMP(I,LMIN) )) THEN          ENDDO
2342           QQ(I,LMIN) = TEMP(I,LMIN)         ENDDO
2343           Q(I,LMIN) = 0.1 * QE(I,LMIN)  c      DO 9220 I = 1,ISTNML
2344          ENDIF  c       IF(IBITSTB(I,LMIN).EQ.1 ) THEN
2345          IF(IBITSTB(I,LMIN).EQ.1 ) P3(I,LMIN) = (2.*B3) *  c        TEMP(I,LMIN) = Q(I,LMIN) / QE(I,LMIN)
2346       1            ( QE(I,LMIN) - Q(I,LMIN) )  c        KH(I,LMIN) = TEMP(I,LMIN) * KH(I,LMIN)
2347   9220  CONTINUE  c        KM(I,LMIN) = TEMP(I,LMIN) * KM(I,LMIN)
2348    c       ENDIF
2349    c       TEMP(I,LMIN) = 0.01 * QQE(I,LMIN)
2350    c       IF((IBITSTB(I,LMIN).EQ.1) .AND.
2351    c    1        ( QQ(I,LMIN) .LE. TEMP(I,LMIN) )) THEN
2352    c        QQ(I,LMIN) = TEMP(I,LMIN)
2353    c        Q(I,LMIN) = 0.1 * QE(I,LMIN)
2354    c       ENDIF
2355    c       IF(IBITSTB(I,LMIN).EQ.1 ) P3(I,LMIN) = (2.*B3) *
2356    c    1            ( QE(I,LMIN) - Q(I,LMIN) )
2357    c9220  CONTINUE
2358           DO L =LMIN,NLEVM1
2359            DO I =1,irun
2360             IF(IBITSTB(I,L).EQ.1 ) THEN
2361              TEMP(I,L) = Q(I,L) / QE(I,L)
2362              KH(I,L) = TEMP(I,L) * KH(I,L)
2363              KM(I,L) = TEMP(I,L) * KM(I,L)
2364             ENDIF
2365             TEMP(I,L) = 0.01 * QQE(I,L)
2366             IF((IBITSTB(I,L).EQ.1) .AND.
2367         1         ( QQ(I,L) .LE. TEMP(I,L) )) THEN
2368              QQ(I,L) = TEMP(I,L)
2369              Q(I,L) = 0.1 * QE(I,L)
2370             ENDIF
2371             IF(IBITSTB(I,L).EQ.1 ) P3(I,L) = (2.*B3) *
2372         1             ( QE(I,L) - Q(I,L) )
2373            ENDDO
2374           ENDDO
2375        ENDIF        ENDIF
2376  C DIMENSIONLESS COEFFS AND P3 (Q GT QE)  C DIMENSIONLESS COEFFS AND P3 (Q GT QE)
2377  C -------------------------------------  C -------------------------------------
# Line 2154  C ------------------------ Line 2386  C ------------------------
2386         ISTNML = irun * ( NLEV - LMIN )         ISTNML = irun * ( NLEV - LMIN )
2387        ENDIF        ENDIF
2388        IF( LMIN .LT. NLEV ) THEN        IF( LMIN .LT. NLEV ) THEN
2389        DO 9122 I =1,ISTNML  c     DO 9122 I =1,ISTNML
2390         P3(I,LMIN) = P3(I,LMIN) * DTAU / XL(I,LMIN)  c      P3(I,LMIN) = P3(I,LMIN) * DTAU / XL(I,LMIN)
2391         TEMP(I,LMIN) = QQE(I,LMIN) * P3(I,LMIN)  c      TEMP(I,LMIN) = QQE(I,LMIN) * P3(I,LMIN)
2392         XQ(I,LMIN) = QQE(I,LMIN) - QQ(I,LMIN)  c      XQ(I,LMIN) = QQE(I,LMIN) - QQ(I,LMIN)
2393  9122  CONTINUE  9122  CONTINUE
2394         DO 9216 I = 1,ISTNML         DO L =LMIN,NLEVM1
2395          IF( ( (IBITSTB(I,LMIN).EQ.1) .AND.          DO I =1,irun
2396       1            ( XQ(I,LMIN) .LT. TEMP(I,LMIN) ) )           P3(I,L) = P3(I,L) * DTAU / XL(I,L)
2397             TEMP(I,L) = QQE(I,L) * P3(I,L)
2398             XQ(I,L) = QQE(I,L) - QQ(I,L)
2399            ENDDO
2400           ENDDO
2401    c      DO 9216 I = 1,ISTNML
2402    c       IF( ( (IBITSTB(I,LMIN).EQ.1) .AND.
2403    c    1            ( XQ(I,LMIN) .LT. TEMP(I,LMIN) ) )
2404    c    2                  .OR.
2405    c    3 ( (IBITSTB(I,LMIN).EQ.0) .AND.
2406    c    4            ( XQ(I,LMIN) .GT. TEMP(I,LMIN) ) ) )
2407    c    5  P3(I,LMIN) = XQ(I,LMIN) / QQE(I,LMIN)
2408    c9216  CONTINUE
2409           DO L =LMIN,NLEVM1
2410            DO I =1,irun
2411             IF( ( (IBITSTB(I,L).EQ.1) .AND.
2412         1             ( XQ(I,L) .LT. TEMP(I,L) ) )
2413       2                  .OR.       2                  .OR.
2414       3 ( (IBITSTB(I,LMIN).EQ.0) .AND.       3       ( (IBITSTB(I,L).EQ.0) .AND.
2415       4            ( XQ(I,LMIN) .GT. TEMP(I,LMIN) ) ) )       4             ( XQ(I,L) .GT. TEMP(I,L) ) ) )
2416       5  P3(I,LMIN) = XQ(I,LMIN) / QQE(I,LMIN)       5   P3(I,L) = XQ(I,L) / QQE(I,L)
2417   9216  CONTINUE          ENDDO
2418           ENDDO
2419        ENDIF        ENDIF
2420   550  CONTINUE   550  CONTINUE
2421  C DIAGNOSTIC PROFILES : INITIAL RI AND QQ  C DIAGNOSTIC PROFILES : INITIAL RI AND QQ
# Line 2181  C -------------------------------------- Line 2430  C --------------------------------------
2430         Z0IN(I) = Z0(I)         Z0IN(I) = Z0(I)
2431         ZETAIN(I) = ZETA(I)         ZETAIN(I) = ZETA(I)
2432  9118  CONTINUE  9118  CONTINUE
2433        DO 9120 I =1,ISTNLV  c     DO 9120 I =1,ISTNLV
2434         RIINIT(I,1) = RI(I,1)  c      RIINIT(I,1) = RI(I,1)
2435         QQINIT(I,1) = QQ(I,1)  c      QQINIT(I,1) = QQ(I,1)
2436  9120  CONTINUE  9120  CONTINUE
2437          DO L =1,NLEV
2438           DO I =1,irun
2439            RIINIT(I,L) = RI(I,L)
2440            QQINIT(I,L) = QQ(I,L)
2441           ENDDO
2442          ENDDO
2443        ENDIF        ENDIF
2444    
2445    C Zero diffusion of TKE above 10 mb
2446          do L = 1,nlev-1
2447            do i = 1,irun
2448    CCC      if(pl(i,L).le.10.) then
2449    CCC       qxlm(i,L) = 0.
2450    CCC      endif
2451             if(pl(i,L).le.150.) then
2452              qxlm(i,L) = min(qxlm(i,L),5.)
2453             endif
2454            enddo
2455          enddo
2456    C
2457  C UPDATE TURBULENT KINETIC ENERGY QQ  C UPDATE TURBULENT KINETIC ENERGY QQ
2458  C ----------------------------------  C ----------------------------------
2459        NLEVMQ = NLEV - LMINQ1        NLEVMQ = NLEV - LMINQ1
2460        ISTNMQ = irun * NLEVMQ        ISTNMQ = irun * NLEVMQ
2461        DO 9306 I =1,ISTNMQ  c     DO 9306 I =1,ISTNMQ
2462        RHOKDZ(I,LMINQ1) = RHODZ1(I,LMINQ1)  c     RHOKDZ(I,LMINQ1) = RHODZ1(I,LMINQ1)
2463       1      * QXLM(I,LMINQ1)  c    1      * QXLM(I,LMINQ1)
2464  9306  CONTINUE  9306  CONTINUE
2465          DO L =LMINQ1,NLEVM1
2466           DO I =1,irun
2467            RHOKDZ(I,L) = RHODZ1(I,L) * QXLM(I,L)
2468           ENDDO
2469          ENDDO
2470         CALL TRBDIF(QQ(1,LMINQ1),P3(1,LMINQ1),RHOKDZ(1,LMINQ1),         CALL TRBDIF(QQ(1,LMINQ1),P3(1,LMINQ1),RHOKDZ(1,LMINQ1),
2471       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)
2472  C  C
2473  C SAVE KH BEFORE ADDING DIMENSIONS FOR USE BY MOIST BOUYANCY CALCULATION  C SAVE KH BEFORE ADDING DIMENSIONS FOR USE BY MOIST BOUYANCY CALCULATION
2474  C  C
2475        DO I = 1,ISTNM1  c     DO I = 1,ISTNM1
2476         KHSAVE(I,1) = KH(I,1)  c      KHSAVE(I,1) = KH(I,1)
2477    c     ENDDO
2478          DO L =1,NLEVM1
2479           DO I =1,irun
2480            KHSAVE(I,L) = KH(I,L)
2481           ENDDO
2482        ENDDO        ENDDO
2483  C  C
2484  C   DIMENSIONAL DIFFUSION COEFFS INCLUDING BACKGROUND AMOUNTS  C   DIMENSIONAL DIFFUSION COEFFS INCLUDING BACKGROUND AMOUNTS
# Line 2208  C Line 2486  C
2486        IF(LMINQ1.GT.1)THEN        IF(LMINQ1.GT.1)THEN
2487         ISTLMQ = irun * (LMINQ1-1)         ISTLMQ = irun * (LMINQ1-1)
2488  CB  CB
2489        DO 9124 I =1,ISTLMQ  c     DO 9124 I =1,ISTLMQ
2490         KM(I,1) = KMBG  c      KM(I,1) = KMBG
2491         KH(I,1) = KHBG  c      KH(I,1) = KHBG
2492  9124  CONTINUE  9124  CONTINUE
2493           DO L =1,LMINQ1-1
2494            DO I =1,irun
2495             KM(I,L) = KMBG
2496             KH(I,L) = KHBG
2497            ENDDO
2498           ENDDO
2499  CE  CE
2500        ENDIF        ENDIF
2501  C  C
2502  CB  CB
2503        DO 9126 I =1,ISTNMQ  c     DO 9126 I =1,ISTNMQ
2504        Q(I,LMINQ1) = 2. * QQ(I,LMINQ1)  c     Q(I,LMINQ1) = 2. * QQ(I,LMINQ1)
2505        Q(I,LMINQ1) = SQRT(Q(I,LMINQ1))  c     Q(I,LMINQ1) = SQRT(Q(I,LMINQ1))
2506        XQ(I,LMINQ1) = XL(I,LMINQ1) * Q(I,LMINQ1)  c     XQ(I,LMINQ1) = XL(I,LMINQ1) * Q(I,LMINQ1)
2507        KM(I,LMINQ1)=XQ(I,LMINQ1)*KM(I,LMINQ1)+KMBG  c     KM(I,LMINQ1)=XQ(I,LMINQ1)*KM(I,LMINQ1)+KMBG
2508        KH(I,LMINQ1)=XQ(I,LMINQ1)*KH(I,LMINQ1)+KHBG  c     KH(I,LMINQ1)=XQ(I,LMINQ1)*KH(I,LMINQ1)+KHBG
2509  9126  CONTINUE  9126  CONTINUE
2510          DO L =LMINQ1,NLEVM1
2511           DO I =1,irun
2512            Q(I,L) = 2. * QQ(I,L)
2513            Q(I,L) = SQRT(Q(I,L))
2514            XQ(I,L)= XL(I,L) * Q(I,L)
2515            KM(I,L)=XQ(I,L)*KM(I,L)+KMBG
2516            KH(I,L)=XQ(I,L)*KH(I,L)+KHBG
2517           ENDDO
2518          ENDDO
2519  CE  CE
2520    C Zero diffusion of u,v,t,q above 10 mb
2521          do L = 1,nlev-1
2522            do i = 1,irun
2523    CCC      if(pl(i,L).le.10.) then
2524    CCC       kh(i,L) = 0.
2525    CCC       km(i,L) = 0.
2526    CCC      endif
2527             if(pl(i,L).le.150.) then
2528              kh(i,L) = min(kh(i,L),5.)
2529              km(i,L) = min(km(i,L),5.)
2530             endif
2531            enddo
2532          enddo
2533  C  C
2534  C   CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: TH AND S  C   CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: TH AND S
2535  C  C
2536        DO 9128 I =1,ISTNLV  c     DO 9128 I =1,ISTNLV
2537        TEMP(I,1) = RHOZPK(I,1) * KH(I,1)  c     TEMP(I,1) = RHOZPK(I,1) * KH(I,1)
2538  9128  CONTINUE  9128  CONTINUE
2539        DO 9130 I =1,ISTNLV  c     DO 9130 I =1,ISTNLV
2540        DELTH(I,1) = 0.  c     DELTH(I,1) = 0.
2541  9130  CONTINUE  9130  CONTINUE
2542          DO L =1,NLEV
2543           DO I =1,irun
2544            TEMP(I,L) = RHOZPK(I,L) * KH(I,L)
2545            DELTH(I,L) = 0.
2546           ENDDO
2547          ENDDO
2548        DO 9132 I =1,irun        DO 9132 I =1,irun
2549        DELTH(I,NLEVP1) = 1.        DELTH(I,NLEVP1) = 1.
2550  9132  CONTINUE  9132  CONTINUE
2551         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)
2552        do i = 1,irun        do i = 1,irun
2553         hsturb(i) = -1.* dths(i)         hsturb(i) = -1.* dths(i)
2554         dhsdtc(i) = -1.* delths(i)         dhsdtc(i) = -1.* delths(i)
# Line 2252  C Line 2564  C
2564         enddo         enddo
2565        enddo        enddo
2566    
2567        DO 9134 I =1,ISTNLV  c     DO 9134 I =1,ISTNLV
2568        RHOKDZ(I,1) = RHODZ2(I,1) * KH(I,1)  c     RHOKDZ(I,1) = RHODZ2(I,1) * KH(I,1)
2569  9134  CONTINUE  9134  CONTINUE
2570        DO 9138 I =1,ISTNLV  c     DO 9138 I =1,ISTNLV
2571        DELSH(I,1) = 0.  c     DELSH(I,1) = 0.
2572  9138  CONTINUE  9138  CONTINUE
2573          DO L =1,NLEV
2574           DO I =1,irun
2575            RHOKDZ(I,L) = RHODZ2(I,L) * KH(I,L)
2576            DELSH(I,L) = 0.
2577           ENDDO
2578          ENDDO
2579        DO 9140 I =1,irun        DO 9140 I =1,irun
2580        DELSH(I,NLEVP1) = 1.        DELSH(I,NLEVP1) = 1.
2581  9140  CONTINUE  9140  CONTINUE
2582    
2583         CALL TRBDIF(SH,DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,NLEV,2,0.,irun)         CALL TRBDIF(SH,DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,NLEV,
2584         .                                     2,0. _d 0,irun)
2585        do i = 1,irun        do i = 1,irun
2586         eturb(i) = -1.* dthl(i)         eturb(i) = -1.* dthl(i)
2587         dedqa(i) = -1.* delthl(i)         dedqa(i) = -1.* delthl(i)
# Line 2282  C Line 2601  C
2601  C   Update Tracers Due to Turbulent Diffusion  C   Update Tracers Due to Turbulent Diffusion
2602  C  C
2603        do i = 1,irun        do i = 1,irun
2604        rhokdz(i,nlev) = 0.0              rhokdz(i,nlev) = 0.0
2605        enddo        enddo
2606    
2607        do nt = 1,ntrace  c     do nt = 1,ntrace
2608        do  i = 1,irun  c     do  i = 1,irun
2609        tracers(i,nlev+1,nt) = tracers(i,nlev,nt)  c     tracers(i,nlev+1,nt) = tracers(i,nlev,nt)
2610        enddo  c     enddo
2611        CALL TRBDIF(tracers(1,1,nt),DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,  c     CALL TRBDIF(tracers(1,1,nt),DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,
2612       .                                                  NLEV,4,0.,irun)  c    .                                         NLEV,4,0. _d 0,irun)
2613        enddo  c     enddo
2614  C  C
2615  C   CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: U AND V  C   CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: U AND V
2616  C  C
2617        DO 9172 I =1,ISTNLV  c     DO 9172 I =1,ISTNLV
2618        RHOKDZ(I,1) = RHODZ2(I,1) * KM(I,1)  c     RHOKDZ(I,1) = RHODZ2(I,1) * KM(I,1)
2619  9172  CONTINUE  9172  CONTINUE
2620         CALL TRBDIF(U,V,RHOKDZ,FLXFAC,DTHS,DELTHS,NLEV,3,0.,irun)        DO L =1,NLEV
2621           DO I =1,irun
2622            RHOKDZ(I,L) = RHODZ2(I,L) * KM(I,L)
2623           ENDDO
2624          ENDDO
2625           CALL TRBDIF(U,V,RHOKDZ,FLXFAC,DTHS,DELTHS,NLEV,3,0. _d 0,irun)
2626  C ( FILL DIAGNOSTIC ARRAYS IF REQUIRED )  C ( FILL DIAGNOSTIC ARRAYS IF REQUIRED )
2627        DO 9174 I =1,ISTNLV  c     DO 9174 I =1,ISTNLV
2628        WU(I,1) = WU(I,1) + RHOKDZ(I,1) * ( U(I,2) - U(I,1) )  c     WU(I,1) = WU(I,1) + RHOKDZ(I,1) * ( U(I,2) - U(I,1) )
2629  9174  CONTINUE  9174  CONTINUE
2630        DO 9176 I =1,ISTNLV  c     DO 9176 I =1,ISTNLV
2631        WV(I,1) = WV(I,1) + RHOKDZ(I,1) * ( V(I,2) - V(I,1) )  c     WV(I,1) = WV(I,1) + RHOKDZ(I,1) * ( V(I,2) - V(I,1) )
2632  9176  CONTINUE  9176  CONTINUE
2633        DO 9300 I = 1,ISTNM1        DO L =1,NLEV
2634        IF ( QQ(I,1) .GT. QQMIN ) THEN         DO I =1,irun
2635           IBITSTB(I,1) = 1          WU(I,L) = WU(I,L) + RHOKDZ(I,L) * ( U(I,L+1) - U(I,L) )
2636           ELSE          WV(I,L) = WV(I,L) + RHOKDZ(I,L) * ( V(I,L+1) - V(I,L) )
2637           IBITSTB(I,1) = 0         ENDDO
2638        ENDIF        ENDDO
2639        IF( IBITSTB(I,1).EQ.1 ) FREQDG(I,1) = FREQDG(I,1) + aitr  c     DO 9300 I = 1,ISTNM1
2640   9300 CONTINUE  c     IF ( QQ(I,1) .GT. QQMIN ) THEN
2641    c        IBITSTB(I,1) = 1
2642    c        ELSE
2643    c        IBITSTB(I,1) = 0
2644    c     ENDIF
2645    c     IF( IBITSTB(I,1).EQ.1 ) FREQDG(I,1) = FREQDG(I,1) + aitr
2646    c9300 CONTINUE
2647          DO L =1,NLEVM1
2648           DO I =1,irun
2649            IF ( QQ(I,L) .GT. QQMIN ) THEN
2650             IBITSTB(I,L) = 1
2651            ELSE
2652             IBITSTB(I,L) = 0
2653            ENDIF
2654            IF( IBITSTB(I,L).EQ.1 ) FREQDG(I,L) = FREQDG(I,L) + aitr
2655           ENDDO
2656          ENDDO
2657         do i = 1,irun         do i = 1,irun
2658          qqcolmin(i) = qq(i,nlev)*0.1          qqcolmin(i) = qq(i,nlev)*0.1
2659          qqcolmax(i) = qq(i,nlev)          qqcolmax(i) = qq(i,nlev)
# Line 2347  C ( FILL DIAGNOSTIC ARRAYS IF REQUIRED ) Line 2687  C ( FILL DIAGNOSTIC ARRAYS IF REQUIRED )
2687        do i=1,irun        do i=1,irun
2688        sz0(i) = sz0(i) + aitr*z0(i)        sz0(i) = sz0(i) + aitr*z0(i)
2689        enddo        enddo
2690        DO I =1,ISTNLV  c     DO I =1,ISTNLV
2691        EU(I,1) = EU(I,1) + AITR*KM(I,1)  c     EU(I,1) = EU(I,1) + AITR*KM(I,1)
2692        enddo  c     enddo
2693        DO I =1,ISTNLV  c     DO I =1,ISTNLV
2694        ET(I,1) = ET(I,1) + AITR*KH(I,1)  c     ET(I,1) = ET(I,1) + AITR*KH(I,1)
2695        enddo  c     enddo
2696          DO L =1,NLEV
2697           DO I =1,irun
2698            EU(I,L) = EU(I,L) + AITR*KM(I,L)
2699            ET(I,L) = ET(I,L) + AITR*KH(I,L)
2700           ENDDO
2701          ENDDO
2702        DO I =1,irun        DO I =1,irun
2703        scu(I) = scu(I) + AITR*cu(I)        scu(I) = scu(I) + AITR*cu(I)
2704        enddo        enddo
# Line 2360  C ( FILL DIAGNOSTIC ARRAYS IF REQUIRED ) Line 2706  C ( FILL DIAGNOSTIC ARRAYS IF REQUIRED )
2706        sct(I) = sct(I) + AITR*ct(I)        sct(I) = sct(I) + AITR*ct(I)
2707        enddo        enddo
2708        IF(tprof) then        IF(tprof) then
2709        do i=1,ISTNM1  c     do i=1,ISTNM1
2710        XLDIAG(I,1) = XLDIAG(I,1) + AITR*XL(I,1)  c     XLDIAG(I,1) = XLDIAG(I,1) + AITR*XL(I,1)
2711        enddo  c     enddo
2712           DO L =1,NLEVM1
2713            DO I =1,irun
2714             XLDIAG(I,L) = XLDIAG(I,L) + AITR*XL(I,L)
2715            ENDDO
2716           ENDDO
2717        endif        endif
2718        FIRST = .FALSE.        FIRST = .FALSE.
2719  C  C
2720  C SAVE XL,CT,XX,YY,ZETA FOR USE BY MOIST BOUYANCY CALCULATION  C SAVE XL,CT,XX,YY,ZETA FOR USE BY MOIST BOUYANCY CALCULATION
2721  C  C
2722        IF(ITER.EQ.ITRTRB)THEN        IF(ITER.EQ.ITRTRB)THEN
2723         DO I = 1,ISTNM1  c      DO I = 1,ISTNM1
2724          XLSAVE(I,1) = XL(I,1)  c       XLSAVE(I,1) = XL(I,1)
2725    c      ENDDO
2726           DO L =1,NLEVM1
2727            DO I =1,irun
2728             XLSAVE(I,L) = XL(I,L)
2729            ENDDO
2730         ENDDO         ENDDO
2731         DO I = 1,irun         DO I = 1,irun
2732          CTSAVE(I) = CT(I)          CTSAVE(I) = CT(I)
# Line 2380  C Line 2736  C
2736         ENDDO         ENDDO
2737        ENDIF        ENDIF
2738    
2739        do i = 1,istnlv  c     do i = 1,istnlv
2740         turbfcc(i,1) = turbfcc(i,1) + fcc(i,1) * aitr  c      turbfcc(i,1) = turbfcc(i,1) + fcc(i,1) * aitr
2741        enddo  c     enddo
2742        do i = 1,irun*nlev  c     do i = 1,irun*nlev
2743         qliq(i,1) = qliq(i,1) + shl(i,1) * aitr  c      qliq(i,1) = qliq(i,1) + shl(i,1) * aitr
2744        enddo  c     enddo
2745          DO L =1,NLEV
2746           DO I =1,irun
2747            turbfcc(I,L) = turbfcc(I,L) + fcc(I,L) * aitr
2748            qliq(I,L) = qliq(I,L) + shl(I,L) * aitr
2749           ENDDO
2750          ENDDO
2751  C  C
2752  C   END OF MAIN LOOP  C   END OF MAIN LOOP
2753  C  C
2754   2000 CONTINUE   2000 CONTINUE
2755        DO 9194 I =1,ISTNLV  c     DO 9194 I =1,ISTNLV
2756         WU(I,1) =  WU(I,1) * AITR  c      WU(I,1) =  WU(I,1) * AITR
2757         WV(I,1) =  WV(I,1) * AITR  c      WV(I,1) =  WV(I,1) * AITR
2758  9194  CONTINUE  9194  CONTINUE
2759          DO L =1,NLEV
2760           DO I =1,irun
2761            WU(I,L) =  WU(I,L) * AITR
2762            WV(I,L) =  WV(I,L) * AITR
2763           ENDDO
2764          ENDDO
2765  C  C
2766        RETURN        RETURN
2767        END        END
# Line 2453  C*************************************** Line 2821  C***************************************
2821    
2822  C Argument List Declarations  C Argument List Declarations
2823        integer nn,n,irun        integer nn,n,irun
2824        real aitr,cp,rgas,undef        _RL aitr,cp,rgas,undef
2825        real VUS(IRUN),VVS(IRUN),VTHV1(IRUN),VTHV2(IRUN)        _RL VUS(IRUN),VVS(IRUN),VTHV1(IRUN),VTHV2(IRUN)
2826        real VTH1(IRUN),VTH2(IRUN),VSH1(IRUN),VSH2(IRUN)        _RL VTH1(IRUN),VTH2(IRUN),VSH1(IRUN),VSH2(IRUN)
2827        real VPK(IRUN),VPKE(IRUN),VPE(IRUN)        _RL VPK(IRUN),VPKE(IRUN),VPE(IRUN)
2828        real VZ0(IRUN),VHS(IRUN),VAHS(IRUN)        _RL VZ0(IRUN),VHS(IRUN),VAHS(IRUN)
2829        integer IVWATER(IRUN)        integer IVWATER(IRUN)
2830        LOGICAL FIRST,LAST        LOGICAL FIRST,LAST
2831        real VRHO(IRUN),VRHOZPK(IRUN)        _RL VRHO(IRUN),VRHOZPK(IRUN)
2832        real VKM(IRUN),VKH(IRUN),VUSTAR(IRUN),VXX(IRUN)        _RL VKM(IRUN),VKH(IRUN),VUSTAR(IRUN),VXX(IRUN)
2833        real VYY(IRUN),VCU(IRUN),VCT(IRUN),VRIB(IRUN)        _RL VYY(IRUN),VCU(IRUN),VCT(IRUN),VRIB(IRUN)
2834        real VZETA(IRUN),VWS(IRUN)        _RL VZETA(IRUN),VWS(IRUN)
2835        real stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)        _RL stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)
2836        real stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)        _RL stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)
2837        LOGICAL LWATER        LOGICAL LWATER
2838        integer IVBITRIB(irun)        integer IVBITRIB(irun)
2839        real VHZ(irun),VPSIM(irun),VAPSIM(irun),VPSIG(irun),VPSIHG(irun)        _RL VHZ(irun),VPSIM(irun),VAPSIM(irun),VPSIG(irun),VPSIHG(irun)
2840        real VTEMP(irun),VDZETA(irun),VDZ0(irun),VDPSIM(irun)        _RL VTEMP(irun),VDZETA(irun),VDZ0(irun),VDPSIM(irun)
2841        real VDPSIH(irun),VZH(irun),VXX0(irun),VYY0(irun)        _RL VDPSIH(irun),VZH(irun),VXX0(irun),VYY0(irun)
2842        real VAPSIHG(irun),VRIB1(irun),VWS1(irun)        _RL VAPSIHG(irun),VRIB1(irun),VWS1(irun)
2843        real VPSIH(irun),VZETAL(irun),VZ0L(irun),VPSIH2(irun),VH0(irun)        _RL VPSIH(irun),VZETAL(irun),VZ0L(irun),VPSIH2(irun),VH0(irun)
2844        real VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)        _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
2845        real VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)        _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
2846        real VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)        _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
2847        real VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)        _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
2848        real VXNUM(irun),VDZETA1(irun),VDZETA2(irun)        _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
2849        real VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)        _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
2850        real VDPSIMC(irun),VDPSIHC(irun)        _RL VDPSIMC(irun),VDPSIHC(irun)
2851    
2852  C Local Variables  C Local Variables
2853        real USTMX3,USTZ0S,Z0MIN,H0BYZ0,USTH0S,H0VEG,Z0VEGM,PRFAC        _RL USTMX3,USTZ0S,Z0MIN,H0BYZ0,USTH0S,H0VEG,Z0VEGM,PRFAC
2854        real XPFAC,DIFSQT        _RL XPFAC,DIFSQT
2855        PARAMETER ( USTMX3 =   0.0632456)        PARAMETER ( USTMX3 =   0.0632456)
2856        PARAMETER ( USTZ0S =   0.2030325E-5)        PARAMETER ( USTZ0S =   0.2030325E-5)
2857        PARAMETER ( Z0MIN  =  USTZ0S/USTMX3)        PARAMETER ( Z0MIN  =  USTZ0S/USTMX3)
2858        PARAMETER ( H0BYZ0 =    30.0    )        PARAMETER ( H0BYZ0 =    30.0    )
2859        PARAMETER ( USTH0S =  H0BYZ0*USTZ0S )        PARAMETER ( USTH0S =  H0BYZ0*USTZ0S )
2860        PARAMETER ( H0VEG  =   0.01     )          PARAMETER ( H0VEG  =   0.01     )
2861        PARAMETER ( Z0VEGM =   0.005    )        PARAMETER ( Z0VEGM =   0.005    )
2862        PARAMETER ( PRFAC  = 0.595864   )        PARAMETER ( PRFAC  = 0.595864   )
2863        PARAMETER ( XPFAC  = .55        )          PARAMETER ( XPFAC  = .55        )
2864        PARAMETER ( DIFSQT  = 3.872983E-3)        PARAMETER ( DIFSQT  = 3.872983E-3)
2865    
2866        real psihdiag(irun),psimdiag(irun)        _RL psihdiag(irun),psimdiag(irun)
2867        real getcon,vk,rvk,vk2,bmdl        _RL getcon,vk,rvk,vk2,bmdl
2868        integer iwater,itype        integer iwater,itype
2869        integer i,iter        integer i,iter
2870  C  C
# Line 2839  C*************************************** Line 3207  C***************************************
3207    
3208  C Argument List Declarations  C Argument List Declarations
3209        integer n,iflag        integer n,iflag
3210        real PHIM(N),PHIH(N),Z(N)        _RL PHIM(N),PHIH(N),Z(N)
3211    
3212  C Local Variables  C Local Variables
3213        integer I1(N),I2(N)        integer I1(N),I2(N)
3214        real ZSTAR(N),E1(N),E2(N),TEMP1(N)        _RL ZSTAR(N),E1(N),E2(N),TEMP1(N)
3215  C  C
3216        real PHIM0(385),ZLINM1(75),ZLINM2(75),ZLINM3(36)        _RL PHIM0(385),ZLINM1(75),ZLINM2(75),ZLINM3(36)
3217        real ZLOGM1(74),ZLOGM2(75),ZLOGM3(50)        _RL ZLOGM1(74),ZLOGM2(75),ZLOGM3(50)
3218        real PHIH0(385),ZLINH1(75),ZLINH2(75),ZLINH3(36)        _RL PHIH0(385),ZLINH1(75),ZLINH2(75),ZLINH3(36)
3219        real ZLOGH1(74),ZLOGH2(75),ZLOGH3(50)        _RL ZLOGH1(74),ZLOGH2(75),ZLOGH3(50)
3220        EQUIVALENCE (PHIM0(1),ZLINM1(1)),(PHIM0(76),ZLINM2(1))        EQUIVALENCE (PHIM0(1),ZLINM1(1)),(PHIM0(76),ZLINM2(1))
3221        EQUIVALENCE (PHIM0(151),ZLINM3(1))        EQUIVALENCE (PHIM0(151),ZLINM3(1))
3222        EQUIVALENCE (PHIM0(187),ZLOGM1(1)),(PHIM0(261),ZLOGM2(1))        EQUIVALENCE (PHIM0(187),ZLOGM1(1)),(PHIM0(261),ZLOGM2(1))
# Line 3139  C*************************************** Line 3507  C***************************************
3507    
3508  C Argument List Declarations  C Argument List Declarations
3509        integer irun,iflag        integer irun,iflag
3510        real VZZ(IRUN),VZH(IRUN),VPSIM(IRUN),VPSIH(IRUN),        _RL VZZ(IRUN),VZH(IRUN),VPSIM(IRUN),VPSIH(IRUN),
3511       1     VX(IRUN),VXS(IRUN),VY(IRUN),VYS(IRUN)       1     VX(IRUN),VXS(IRUN),VY(IRUN),VYS(IRUN)
3512    
3513  C Local Variables  C Local Variables
3514        real ZWM,RZWM,Z0M,ZCM,RZCM,CM1,CM2,CM6,CM7,CM8ARG,YCM        _RL ZWM,RZWM,Z0M,ZCM,RZCM,CM1,CM2,CM6,CM7,CM8ARG,YCM
3515        PARAMETER ( ZWM     =    1.    )        PARAMETER ( ZWM     =    1.    )
3516        PARAMETER ( RZWM    =  1./ZWM  )        PARAMETER ( RZWM    =  1./ZWM  )
3517        PARAMETER ( Z0M     =    0.2    )        PARAMETER ( Z0M     =    0.2    )
# Line 3157  C Local Variables Line 3525  C Local Variables
3525        PARAMETER ( YCM     =  6. / ( 1. + 6.*CM1*ZCM )  )        PARAMETER ( YCM     =  6. / ( 1. + 6.*CM1*ZCM )  )
3526    
3527        integer INTSTB(irun),INTZ0(irun)        integer INTSTB(irun),INTZ0(irun)
3528        real ZZ0(irun),Z(irun),Z2(irun),Z1(irun),Z0(irun)        _RL ZZ0(irun),Z(irun),Z2(irun),Z1(irun),Z0(irun)
3529        real X0(irun),X1(irun),Y0(irun),Y1(irun)        _RL X0(irun),X1(irun),Y0(irun),Y1(irun)
3530        real PSI2(irun),TEMP(irun)        _RL PSI2(irun),TEMP(irun)
3531        real HZ(irun),ARG0(irun),ARG1(irun),DX(irun)        _RL HZ(irun),ARG0(irun),ARG1(irun),DX(irun)
3532        real X0NUM(irun),X1NUM(irun),X0DEN(irun)        _RL X0NUM(irun),X1NUM(irun),X0DEN(irun)
3533        real X1DEN(irun),Y1DEN(irun),Z2ZWM(irun)        _RL X1DEN(irun),Y1DEN(irun),Z2ZWM(irun)
3534        real cm3,cm4,cm5,cm8        _RL cm3,cm4,cm5,cm8
3535        integer ibit,index        integer ibit,indx
3536        integer i        integer i
3537  C  C
3538        CM3 =   sqrt( 0.2/CM1-0.01 )        CM3 =   sqrt( 0.2/CM1-0.01 )
# Line 3199  C ************************************** Line 3567  C **************************************
3567  C  C
3568        IF(IBIT.LE.0)  GO TO 100        IF(IBIT.LE.0)  GO TO 100
3569  C  C
3570        INDEX = 0        indx = 0
3571        DO 9002 I = 1,IRUN        DO 9002 I = 1,IRUN
3572         IF (INTSTB(I).EQ.1)THEN         IF (INTSTB(I).EQ.1)THEN
3573          INDEX = INDEX + 1          indx = indx + 1
3574          Z(INDEX) = VZZ(I)          Z(indx) = VZZ(I)
3575          Z0(INDEX) = ZZ0(I)          Z0(indx) = ZZ0(I)
3576         ENDIF         ENDIF
3577   9002 CONTINUE   9002 CONTINUE
3578  C  C
# Line 3212  C Line 3580  C
3580         Z(I) = -18. * Z(I)         Z(I) = -18. * Z(I)
3581         Z0(I) = -18. * Z0(I)         Z0(I) = -18. * Z0(I)
3582   9004 CONTINUE   9004 CONTINUE
3583    
3584        CALL PHI( Z,X1,Y1,IFLAG,IBIT )        CALL PHI( Z,X1,Y1,IFLAG,IBIT )
3585        CALL PHI( Z0,X0,Y0,IFLAG,IBIT )        CALL PHI( Z0,X0,Y0,IFLAG,IBIT )
3586    
3587  C ****************************  C ****************************
3588  C *****    COMPUTE PSIM  *****  C *****    COMPUTE PSIM  *****
3589  C ****************************  C ****************************
# Line 3242  C Line 3610  C
3610         PSI2(I) = PSI2(I) + DX(I)         PSI2(I) = PSI2(I) + DX(I)
3611   9006 CONTINUE   9006 CONTINUE
3612  C  C
3613        INDEX = 0        indx = 0
3614        DO 9008 I = 1,IRUN        DO 9008 I = 1,IRUN
3615         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3616          INDEX = INDEX + 1          indx = indx + 1
3617          VPSIM(I) = PSI2(INDEX)          VPSIM(I) = PSI2(indx)
3618          VX(I) = X1(INDEX)          VX(I) = X1(indx)
3619          VXS(I) = X0(INDEX)          VXS(I) = X0(indx)
3620         ENDIF         ENDIF
3621   9008 CONTINUE   9008 CONTINUE
3622  C  C
# Line 3275  C Line 3643  C
3643         PSI2(I) = PSI2(I) - Y1(I) + Y0(I)         PSI2(I) = PSI2(I) - Y1(I) + Y0(I)
3644   9010 CONTINUE   9010 CONTINUE
3645  C  C
3646        INDEX = 0        indx = 0
3647        DO 9012 I = 1,IRUN        DO 9012 I = 1,IRUN
3648         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3649         INDEX = INDEX + 1         indx = indx + 1
3650         VPSIH(I) = PSI2(INDEX)         VPSIH(I) = PSI2(indx)
3651         VY(I) = Y1(INDEX)         VY(I) = Y1(indx)
3652         VYS(I) = Y0(INDEX)         VYS(I) = Y0(indx)
3653         ENDIF         ENDIF
3654   9012 CONTINUE   9012 CONTINUE
3655  C  C
# Line 3304  C Line 3672  C
3672         ENDIF         ENDIF
3673   9014 CONTINUE   9014 CONTINUE
3674        IF(IBIT.LE.0)  GO TO 300        IF(IBIT.LE.0)  GO TO 300
3675        INDEX = 0        indx = 0
3676  #ifdef CRAY  #ifdef CRAY
3677  CDIR$ NOVECTOR  CDIR$ NOVECTOR
3678  #endif  #endif
3679        DO 9016 I = 1,IRUN        DO 9016 I = 1,IRUN
3680         IF (INTSTB(I).EQ.1)THEN         IF (INTSTB(I).EQ.1)THEN
3681          INDEX = INDEX + 1          indx = indx + 1
3682          Z(INDEX) = VZZ(I)          Z(indx) = VZZ(I)
3683          Z0(INDEX) = ZZ0(I)          Z0(indx) = ZZ0(I)
3684          ARG1(INDEX) = VZH(I)          ARG1(indx) = VZH(I)
3685         ENDIF         ENDIF
3686   9016 CONTINUE   9016 CONTINUE
3687  #ifdef CRAY  #ifdef CRAY
# Line 3370  C Line 3738  C
3738         PSI2(I) = TEMP(I) + CM6 * ARG1(I)         PSI2(I) = TEMP(I) + CM6 * ARG1(I)
3739   9020 CONTINUE   9020 CONTINUE
3740  C  C
3741        INDEX = 0        indx = 0
3742        DO 9030 I = 1,IRUN        DO 9030 I = 1,IRUN
3743         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3744         INDEX = INDEX + 1         indx = indx + 1
3745         VPSIM(I) = PSI2(INDEX)         VPSIM(I) = PSI2(indx)
3746         VX(I) = X1(INDEX)         VX(I) = X1(indx)
3747         VXS(I) = X0(INDEX)         VXS(I) = X0(indx)
3748         ENDIF         ENDIF
3749   9030 CONTINUE   9030 CONTINUE
3750  C  C
# Line 3402  C Line 3770  C
3770         PSI2(I) = TEMP(I) + ARG0(I) * ARG1(I)         PSI2(I) = TEMP(I) + ARG0(I) * ARG1(I)
3771   9024 CONTINUE   9024 CONTINUE
3772  C  C
3773        INDEX = 0        indx = 0
3774        DO 9026 I = 1,IRUN        DO 9026 I = 1,IRUN
3775         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3776         INDEX = INDEX + 1         indx = indx + 1
3777         VPSIH(I) = PSI2(INDEX)         VPSIH(I) = PSI2(indx)
3778         VY(I) = Y1(INDEX)         VY(I) = Y1(indx)
3779         VYS(I) = X0(INDEX)         VYS(I) = X0(indx)
3780         ENDIF         ENDIF
3781   9026 CONTINUE   9026 CONTINUE
3782  C  C
# Line 3458  C*************************************** Line 3826  C***************************************
3826    
3827  C Argument List Declarations  C Argument List Declarations
3828        integer irun,nlev,init,lmin,lminq,lminq1        integer irun,nlev,init,lmin,lminq,lminq1
3829        real cp        _RL cp
3830        real STRT(irun,NLEV),DW2(irun,NLEV),DZ3(irun,NLEV)        _RL STRT(irun,NLEV),DW2(irun,NLEV),DZ3(irun,NLEV)
3831        real Q(irun,NLEV),VKZM(irun,NLEV-1),VKZE(irun,NLEV-1)        _RL Q(irun,NLEV),VKZM(irun,NLEV-1),VKZE(irun,NLEV-1)
3832        real DTHV(irun,NLEV),DPK(irun,NLEV),DU(irun,NLEV)        _RL DTHV(irun,NLEV),DPK(irun,NLEV),DU(irun,NLEV)
3833        real DV(irun,NLEV)        _RL DV(irun,NLEV)
3834        real QXLM(irun,NLEV-1),XL(irun,NLEV-1)        _RL QXLM(irun,NLEV-1),XL(irun,NLEV-1)
3835        real DZITRP(irun,nlev-1),STBFCN(irun,nlev)        _RL DZITRP(irun,nlev-1),STBFCN(irun,nlev)
3836        real XL0(irun,nlev),Q1(irun,nlev-1)        _RL XL0(irun,nlev),Q1(irun,nlev-1)
3837        real WRKIT1(irun,nlev-1),WRKIT2(irun,nlev-1)        _RL WRKIT1(irun,nlev-1),WRKIT2(irun,nlev-1)
3838        real WRKIT3(irun,nlev-1)        _RL WRKIT3(irun,nlev-1)
3839        real WRKIT4(irun,nlev-1)        _RL WRKIT4(irun,nlev-1)
3840        INTEGER INT1(irun,nlev), INT2(irun,nlev-1)        INTEGER INT1(irun,nlev), INT2(irun,nlev-1)
3841    
3842  C Local Variables  C Local Variables
3843        real rf1,rf2,e5,d4,d1,rfc,ricr,alpha,dzcnv,xl0cnv,xl0min        _RL rf1,rf2,e5,d4,d1,rfc,ricr,alpha,dzcnv,xl0cnv,xl0min
3844        real clmt,clmt53        _RL clmt,clmt53
3845        PARAMETER ( RF1     = 0.2340678 )        PARAMETER ( RF1     = 0.2340678 )
3846        PARAMETER ( RF2     = 0.2231172 )        PARAMETER ( RF2     = 0.2231172 )
3847        PARAMETER ( E5      = 49.66     )          PARAMETER ( E5      = 49.66     )
3848        PARAMETER ( D4      = 2.6532122E-2 )        PARAMETER ( D4      = 2.6532122E-2 )
3849        PARAMETER ( D1      = D4 * E5   )        PARAMETER ( D1      = D4 * E5   )
3850        PARAMETER ( RFC     = 0.1912323 )        PARAMETER ( RFC     = 0.1912323 )
# Line 3485  C Local Variables Line 3853  C Local Variables
3853        PARAMETER ( DZCNV   = 100.      )        PARAMETER ( DZCNV   = 100.      )
3854        PARAMETER ( XL0CNV  = DZCNV * ALPHA )        PARAMETER ( XL0CNV  = DZCNV * ALPHA )
3855        PARAMETER ( XL0MIN  = 1.        )        PARAMETER ( XL0MIN  = 1.        )
3856        PARAMETER ( CLMT    = 0.23      )          PARAMETER ( CLMT    = 0.23      )
3857        PARAMETER ( CLMT53  = 5. * CLMT / 3. )        PARAMETER ( CLMT53  = 5. * CLMT / 3. )
3858    
3859        integer ibit,nlevm1,nlevp1,istnlv,istnm1,nlevml,istnml,Lp1        integer ibit,nlevm1,nlevp1,istnlv,istnm1,nlevml,istnml,Lp1
# Line 3501  C Line 3869  C
3869  C  C
3870  C COMPUTE DEPTHS OF UNSTABLE LAYERS  C COMPUTE DEPTHS OF UNSTABLE LAYERS
3871  C =================================  C =================================
3872        DO 10 I=1,ISTNLV  c     DO 10 I=1,ISTNLV
3873            STBFCN(I,1) = STRT(I,1) - RICR * DW2(I,1)  c         STBFCN(I,1) = STRT(I,1) - RICR * DW2(I,1)
3874                                  INT1(I,1) = 0  c                               INT1(I,1) = 0
3875        IF( STBFCN(I,1).LE.0. ) INT1(I,1) = 1  c     IF( STBFCN(I,1).LE.0. ) INT1(I,1) = 1
3876     10 CONTINUE  c  10 CONTINUE
3877        DO 20 I=1,ISTNM1        DO L =1,NLEV
3878                                                  INT2(I,1) = 0         DO I =1,irun
3879        IF( INT1(I,1).EQ.1 .XOR. INT1(I,2).EQ.1 ) INT2(I,1) = 1           STBFCN(I,L) = STRT(I,L) - RICR * DW2(I,L)
3880     20 CONTINUE           INT1(I,L) = 0
3881             IF( STBFCN(I,L).LE.0. ) INT1(I,L) = 1
3882           ENDDO
3883          ENDDO
3884    
3885    c     DO 20 I=1,ISTNM1
3886    c                                               INT2(I,1) = 0
3887    c     IF( (INT1(I,1).EQ.1) .NEQV. (INT1(I,2).EQ.1) ) INT2(I,1) = 1
3888    c  20 CONTINUE
3889          DO L =1,NLEVM1
3890           DO I =1,irun
3891            INT2(I,L) = 0
3892            IF( (INT1(I,L).EQ.1) .NEQV. (INT1(I,L+1).EQ.1) ) INT2(I,L) = 1
3893           ENDDO
3894          ENDDO
3895  C  C
3896        DO 40 LMIN = 1,NLEV        DO 40 LMIN = 1,NLEV
3897        IBIT = 0        IBIT = 0
# Line 3522  C Line 3904  C
3904     50 CONTINUE     50 CONTINUE
3905        LMIN = 1        LMIN = 1
3906  C  C
3907        DO 60 I=1,ISTNM1  c     DO 60 I=1,ISTNM1
3908        XL0(I,1) = 0.  c     XL0(I,1) = 0.
3909     60 CONTINUE  c  60 CONTINUE
3910          DO L =1,NLEVM1
3911           DO I =1,irun
3912            XL0(I,L) = 0.
3913           ENDDO
3914          ENDDO
3915        DO 70 I=1,irun        DO 70 I=1,irun
3916        XL0(I,NLEV) = DZ3(I,NLEV)        XL0(I,NLEV) = DZ3(I,NLEV)
3917     70 CONTINUE     70 CONTINUE
# Line 3540  C Line 3927  C
3927       .              WRKIT1,WRKIT2,WRKIT3,WRKIT4,CP,irun )       .              WRKIT1,WRKIT2,WRKIT3,WRKIT4,CP,irun )
3928        LP1 = LMIN1 + 1        LP1 = LMIN1 + 1
3929  C  C
3930        DO 80 I=1,ISTNML  c     DO 80 I=1,ISTNML
3931                                                       INT2(I,LMIN1) = 0  c                                                    INT2(I,LMIN1) = 0
3932        IF( INT1(I,LMIN1).EQ.1 .OR. INT1(I,LP1).EQ.1 ) INT2(I,LMIN1) = 1  c     IF( INT1(I,LMIN1).EQ.1 .OR. INT1(I,LP1).EQ.1 ) INT2(I,LMIN1) = 1
3933        IF( INT2(I,LMIN1).EQ.1 )  c     IF( INT2(I,LMIN1).EQ.1 )
3934       .     XL0(I,LMIN1) = (0.5+DZITRP(I,LMIN1)) * DZ3(I,LP1)  c    .     XL0(I,LMIN1) = (0.5+DZITRP(I,LMIN1)) * DZ3(I,LP1)
3935     80 CONTINUE  c  80 CONTINUE
3936          DO L =LMIN1,NLEVM1
3937           DO I =1,irun
3938            INT2(I,L) = 0
3939            IF ( INT1(I,L).EQ.1 .OR. INT1(I,L+1).EQ.1 ) INT2(I,L) = 1
3940            IF ( INT2(I,L).EQ.1 )
3941         .       XL0(I,L) = (0.5+DZITRP(I,L)) * DZ3(I,L+1)
3942           ENDDO
3943          ENDDO
3944        DO 90 I=1,irun        DO 90 I=1,irun
3945        INT2(I,NLEVM1) = INT1(I,NLEV)        INT2(I,NLEVM1) = INT1(I,NLEV)
3946     90 CONTINUE     90 CONTINUE
3947  C  C
3948        DO 100 I=1,ISTNML  c     DO 100 I=1,ISTNML
3949        IF( INT2(I,LMIN1).EQ.1 ) THEN  c     IF( INT2(I,LMIN1).EQ.1 ) THEN
3950        XL0(I,LP1) = XL0(I,LP1) + ( (0.5-DZITRP(I,LMIN1)) * DZ3(I,LP1) )  c     XL0(I,LP1) = XL0(I,LP1) + ( (0.5-DZITRP(I,LMIN1)) * DZ3(I,LP1) )
3951        ENDIF  c     ENDIF
3952    100 CONTINUE  c 100 CONTINUE
3953          DO L =LMIN1,NLEVM1
3954           DO I =1,irun
3955            IF ( INT2(I,L).EQ.1 ) THEN
3956             XL0(I,L+1) = XL0(I,L+1) + ( (0.5-DZITRP(I,L)) * DZ3(I,L+1) )
3957            ENDIF
3958           ENDDO
3959          ENDDO
3960        IF (LMIN.GT.1) GOTO 400        IF (LMIN.GT.1) GOTO 400
3961        DO 110 I=1,irun        DO 110 I=1,irun
3962        IF( INT1(I,1).EQ.1 ) XL0(I,1) = XL0(I,1) + DZ3(I,1)        IF( INT1(I,1).EQ.1 ) XL0(I,1) = XL0(I,1) + DZ3(I,1)
# Line 3593  C Line 3995  C
3995   1000 CONTINUE   1000 CONTINUE
3996   1100 CONTINUE   1100 CONTINUE
3997  C  C
3998        DO 150 I = 1,ISTNLV  c     DO 150 I = 1,ISTNLV
3999        IF( XL0(I,1).LT.XL0CNV ) XL0(I,1) = XL0CNV  c     IF( XL0(I,1).LT.XL0CNV ) XL0(I,1) = XL0CNV
4000    150 CONTINUE  c 150 CONTINUE
4001          DO L =1,NLEV
4002           DO I =1,irun
4003            IF ( XL0(I,L).LT.XL0CNV ) XL0(I,L) = XL0CNV
4004           ENDDO
4005          ENDDO
4006  C  C
4007  C *********************************************************************  C *********************************************************************
4008  C ****          DETERMINE MIXING LENGTHS FOR STABLE LAYERS          ***  C ****          DETERMINE MIXING LENGTHS FOR STABLE LAYERS          ***
# Line 3605  C Line 4012  C
4012  C  C
4013        IF(LMINQ.GT.1) THEN        IF(LMINQ.GT.1) THEN
4014                      ISTLMQ = irun * LMINQ1                      ISTLMQ = irun * LMINQ1
4015        DO 160 I  = 1,ISTLMQ  c     DO 160 I  = 1,ISTLMQ
4016        INT2(I,1) = 1 - INT1(I,1)  c     INT2(I,1) = 1 - INT1(I,1)
4017    160 CONTINUE  c 160 CONTINUE
4018          DO L =1,LMINQ1
4019           DO I =1,irun
4020            INT2(I,L) = 1 - INT1(I,L)
4021           ENDDO
4022          ENDDO
4023        ENDIF        ENDIF
4024        IF(LMINQ.LT.NLEV) THEN        IF(LMINQ.LT.NLEV) THEN
4025          ISTNMQ = irun * (NLEV-LMINQ)          ISTNMQ = irun * (NLEV-LMINQ)
4026        DO 170 I = 1,ISTNMQ  c     DO 170 I = 1,ISTNMQ
4027        IF( INT1(I,LMINQ).EQ.0 ) THEN  c     IF( INT1(I,LMINQ).EQ.0 ) THEN
4028             XL0(I,LMINQ) =         Q(I,LMINQ) / XL0(I,LMINQ)  c          XL0(I,LMINQ) =         Q(I,LMINQ) / XL0(I,LMINQ)
4029             XL0(I,LMINQ) =       XL0(I,LMINQ) * XL0(I,LMINQ) + 1.0E-20  c          XL0(I,LMINQ) =       XL0(I,LMINQ) * XL0(I,LMINQ) + 1.0E-20
4030             XL0(I,LMINQ) =    STBFCN(I,LMINQ) + XL0(I,LMINQ)  c          XL0(I,LMINQ) =    STBFCN(I,LMINQ) + XL0(I,LMINQ)
4031             XL0(I,LMINQ) = SQRT( XL0(I,LMINQ) )  c          XL0(I,LMINQ) = SQRT( XL0(I,LMINQ) )
4032             XL0(I,LMINQ) =         Q(I,LMINQ) / XL0(I,LMINQ)  c          XL0(I,LMINQ) =         Q(I,LMINQ) / XL0(I,LMINQ)
4033        ENDIF  c     ENDIF
4034                                     INT2(I,LMINQ)  = 0  c                                  INT2(I,LMINQ)  = 0
4035        IF( XL0(I,LMINQ).LT.XL0MIN ) INT2(I,LMINQ)  = 1  c     IF( XL0(I,LMINQ).LT.XL0MIN ) INT2(I,LMINQ)  = 1
4036    170 CONTINUE  c 170 CONTINUE
4037           DO L =LMINQ,NLEVM1
4038            DO I =1,irun
4039             IF( INT1(I,L).EQ.0 ) THEN
4040               XL0(I,L) =         Q(I,L) / XL0(I,L)
4041               XL0(I,L) =       XL0(I,L) * XL0(I,L) + 1.0E-20
4042               XL0(I,L) =    STBFCN(I,L) + XL0(I,L)
4043               XL0(I,L) = SQRT( XL0(I,L) )
4044               XL0(I,L) =         Q(I,L) / XL0(I,L)
4045             ENDIF
4046             INT2(I,L)  = 0
4047             IF( XL0(I,L).LT.XL0MIN ) INT2(I,L)  = 1
4048            ENDDO
4049           ENDDO
4050        ENDIF        ENDIF
4051  C  C
4052   1200 CONTINUE   1200 CONTINUE
4053  C  C
4054        IF(INIT.EQ.2) THEN        IF(INIT.EQ.2) THEN
4055        DO 180 I = 1,ISTNM1  c     DO 180 I = 1,ISTNM1
4056        INT2(I,1) = 1 - INT1(I,1)  c     INT2(I,1) = 1 - INT1(I,1)
4057    180 CONTINUE  c 180 CONTINUE
4058        ENDIF         DO L =1,NLEVM1
4059        DO 190 I = 1,ISTNM1          DO I =1,irun
4060        IF( INT2(I,1).EQ.1 ) XL0(I,1) = XL0MIN           INT2(I,L) = 1 - INT1(I,L)
4061    190 CONTINUE           IF ( INT2(I,L).EQ.1 ) XL0(I,L) = XL0MIN
4062            ENDDO
4063           ENDDO
4064          ENDIF
4065    c     DO 190 I = 1,ISTNM1
4066    c     IF( INT2(I,1).EQ.1 ) XL0(I,1) = XL0MIN
4067    c 190 CONTINUE
4068          DO L =1,NLEVM1
4069           DO I =1,irun
4070            IF ( INT2(I,L).EQ.1 ) XL0(I,L) = XL0MIN
4071           ENDDO
4072          ENDDO
4073  C  C
4074  C *********************************************************************  C *********************************************************************
4075  C ****             LENGTH SCALE XL FROM XL0 AND VKZE              ****  C ****             LENGTH SCALE XL FROM XL0 AND VKZE              ****
# Line 3641  C ************************************** Line 4077  C **************************************
4077  C  C
4078   1400 CONTINUE   1400 CONTINUE
4079  C  C
4080        DO 200 I = 1,ISTNM1  c     DO 200 I = 1,ISTNM1
4081        XL(I,1) = XL0(I,1) * VKZE(I,1) / ( XL0(I,1)+VKZE(I,1) )  c     XL(I,1) = XL0(I,1) * VKZE(I,1) / ( XL0(I,1)+VKZE(I,1) )
4082    200 CONTINUE  c 200 CONTINUE
4083          DO L =1,NLEVM1
4084           DO I =1,irun
4085            XL(I,L) = XL0(I,L) * VKZE(I,L) / ( XL0(I,L)+VKZE(I,L) )
4086           ENDDO
4087          ENDDO
4088  C  C
4089  C *********************************************************************  C *********************************************************************
4090  C ****       CLMT53 TIMES Q TIMES LENGTH SCALE AT MID LEVELS        ***  C ****       CLMT53 TIMES Q TIMES LENGTH SCALE AT MID LEVELS        ***
# Line 3651  C ************************************** Line 4092  C **************************************
4092  C  C
4093        IF(INIT.EQ.1) GOTO 1700        IF(INIT.EQ.1) GOTO 1700
4094                     ISTNMQ = irun * (NLEV-LMINQ1)                     ISTNMQ = irun * (NLEV-LMINQ1)
4095        DO 210 I = 1,ISTNMQ  c     DO 210 I = 1,ISTNMQ
4096                                                Q1(I,LMINQ1) = Q(I,LMINQ1)  c                                             Q1(I,LMINQ1) = Q(I,LMINQ1)
4097                                              INT1(I,LMINQ1) = 0  c                                           INT1(I,LMINQ1) = 0
4098        IF(    Q(I,LMINQ1).LE.Q(I,LMINQ1+1) ) INT1(I,LMINQ1) = 1  c     IF(    Q(I,LMINQ1).LE.Q(I,LMINQ1+1) ) INT1(I,LMINQ1) = 1
4099        IF( INT1(I,LMINQ1).EQ.1 )  THEN  c     IF( INT1(I,LMINQ1).EQ.1 )  THEN
4100             XL0(I,LMINQ1) = XL0(I,LMINQ1+1)  c          XL0(I,LMINQ1) = XL0(I,LMINQ1+1)
4101              Q1(I,LMINQ1) =   Q(I,LMINQ1+1)  c           Q1(I,LMINQ1) =   Q(I,LMINQ1+1)
4102        ENDIF  c     ENDIF
4103    210 CONTINUE  c 210 CONTINUE
4104  C        DO L =LMINQ1,NLEVM1
4105        DO 240 I = 1,ISTNMQ         DO I =1,irun
4106        QXLM(I,LMINQ1) =   XL0(I,LMINQ1)*VKZM(I,LMINQ1)          Q1(I,L) = Q(I,L)
4107       .               / ( XL0(I,LMINQ1)+VKZM(I,LMINQ1) )          INT1(I,L) = 0
4108        QXLM(I,LMINQ1) = CLMT53 * Q1(I,LMINQ1)*QXLM(I,LMINQ1)          IF (    Q(I,L).LE.Q(I,L+1) ) INT1(I,L) = 1
4109    240 CONTINUE          IF ( INT1(I,L).EQ.1 )  THEN
4110               XL0(I,L) = XL0(I,L+1)
4111                Q1(I,L) =   Q(I,L+1)
4112            ENDIF
4113           ENDDO
4114          ENDDO
4115    C
4116    c     DO 240 I = 1,ISTNMQ
4117    c     QXLM(I,LMINQ1) =   XL0(I,LMINQ1)*VKZM(I,LMINQ1)
4118    c    .               / ( XL0(I,LMINQ1)+VKZM(I,LMINQ1) )
4119    c     QXLM(I,LMINQ1) = CLMT53 * Q1(I,LMINQ1)*QXLM(I,LMINQ1)
4120    c 240 CONTINUE
4121          DO L =LMINQ1,NLEVM1
4122           DO I =1,irun
4123            QXLM(I,L) =   XL0(I,L)*VKZM(I,L)
4124         .                 / ( XL0(I,L)+VKZM(I,L) )
4125            QXLM(I,L) = CLMT53 * Q1(I,L)*QXLM(I,L)
4126           ENDDO
4127          ENDDO
4128  C  C
4129   1700 CONTINUE   1700 CONTINUE
4130  C  C
# Line 3698  C*************************************** Line 4157  C***************************************
4157    
4158  C Argument List Declarations  C Argument List Declarations
4159        integer irun,nlev        integer irun,nlev
4160        real cp        _RL cp
4161        real STBFCN(irun,NLEV+1)        _RL STBFCN(irun,NLEV+1)
4162        integer INTCHG(irun,NLEV)        integer INTCHG(irun,NLEV)
4163        real DTHV(irun,NLEV+1),DPK(irun,NLEV+1)        _RL DTHV(irun,NLEV+1),DPK(irun,NLEV+1)
4164        real DU(irun,NLEV+1),DV(irun,NLEV+1)        _RL DU(irun,NLEV+1),DV(irun,NLEV+1)
4165        real DZITRP(irun,NLEV+1)        _RL DZITRP(irun,NLEV+1)
4166        real AAA(irun,NLEV),BBB(irun,NLEV)        _RL AAA(irun,NLEV),BBB(irun,NLEV)
4167        real CCC(irun,NLEV),DDD(irun,NLEV)        _RL CCC(irun,NLEV),DDD(irun,NLEV)
4168    
4169  C Local Variables  C Local Variables
4170        real rf1,rf2,e5,d4,d1,rfc,ricr        _RL rf1,rf2,e5,d4,d1,rfc,ricr
4171        PARAMETER ( RF1     = 0.2340678 )        PARAMETER ( RF1     = 0.2340678 )
4172        PARAMETER ( RF2     = 0.2231172 )        PARAMETER ( RF2     = 0.2231172 )
4173        PARAMETER ( E5      = 49.66     )          PARAMETER ( E5      = 49.66     )
4174        PARAMETER ( D4      = 2.6532122E-2 )        PARAMETER ( D4      = 2.6532122E-2 )
4175        PARAMETER ( D1      = D4 * E5   )        PARAMETER ( D1      = D4 * E5   )
4176        PARAMETER ( RFC     = 0.1912323 )        PARAMETER ( RFC     = 0.1912323 )
4177        PARAMETER ( RICR    = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 )  )        PARAMETER ( RICR    = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 )  )
4178    
4179        integer istnlv        integer istnlv
4180        integer i        integer I,L
4181  C  C
4182  C *********************************************************************  C *********************************************************************
4183  C ****           QUADRATIC INTERPOLATION OF RI TO RICR VIA          ***  C ****           QUADRATIC INTERPOLATION OF RI TO RICR VIA          ***
# Line 3726  C ****           LINEAR INTERPOLATION OF Line 4185  C ****           LINEAR INTERPOLATION OF
4185  C *********************************************************************  C *********************************************************************
4186  C  C
4187        ISTNLV = irun*NLEV        ISTNLV = irun*NLEV
4188        DO 10 I=1,ISTNLV  c     DO 10 I=1,ISTNLV
4189        DZITRP(I,1) = 0.  c     DZITRP(I,1) = 0.
4190     10 CONTINUE  c  10 CONTINUE
4191        DO 20 I=1,ISTNLV        DO L =1,NLEV
4192        IF( INTCHG(I,1).EQ.1 ) THEN         DO I =1,irun
4193        DDD(I,1) = (         CP  *(DTHV(I,2)*DPK(I,1)          DZITRP(I,L) = 0.
4194       .                         + DTHV(I,1)*DPK(I,2)) )         ENDDO
4195       .         - ( (2.*RICR)  * ( DU(I,2)* DU(I,1)        ENDDO
4196       .                         +   DV(I,2)* DV(I,1)) )  c     DO 20 I=1,ISTNLV
4197        AAA(I,1) = STBFCN(I,1) + STBFCN(I,2)  c     IF( INTCHG(I,1).EQ.1 ) THEN
4198        BBB(I,1) = STBFCN(I,1) - STBFCN(I,2)  c     DDD(I,1) = (         CP  *(DTHV(I,2)*DPK(I,1)
4199        CCC(I,1) =            1. / BBB(I,1)  c    .                         + DTHV(I,1)*DPK(I,2)) )
4200        DZITRP(I,1) = AAA(I,1) * CCC(I,1)  c    .         - ( (2.*RICR)  * ( DU(I,2)* DU(I,1)
4201           AAA(I,1) = AAA(I,1) - DDD(I,1)  c    .                         +   DV(I,2)* DV(I,1)) )
4202           DDD(I,1) =        (   DDD(I,1) *    DDD(I,1) )  c     AAA(I,1) = STBFCN(I,1) + STBFCN(I,2)
4203       .            - 4. * (STBFCN(I,2) * STBFCN(I,1) )  c     BBB(I,1) = STBFCN(I,1) - STBFCN(I,2)
4204           DDD(I,1) = DDD(I,1)*CCC(I,1)*CCC(I,1)  c     CCC(I,1) =            1. / BBB(I,1)
4205           DDD(I,1) =    SQRT( DDD(I,1) )  c     DZITRP(I,1) = AAA(I,1) * CCC(I,1)
4206        ENDIF  c        AAA(I,1) = AAA(I,1) - DDD(I,1)
4207    c        DDD(I,1) =        (   DDD(I,1) *    DDD(I,1) )
4208    c    .            - 4. * (STBFCN(I,2) * STBFCN(I,1) )
4209    c        DDD(I,1) = DDD(I,1)*CCC(I,1)*CCC(I,1)
4210    c        DDD(I,1) =    SQRT( DDD(I,1) )
4211    c     ENDIF
4212    C
4213    c     IF( INTCHG(I,1).EQ.1 .AND. AAA(I,1).NE.0. ) THEN
4214    c     DZITRP(I,1) = ( BBB(I,1)*(1.-DDD(I,1)) ) / AAA(I,1)
4215    c     ENDIF
4216    C
4217    c     DZITRP(I,1) = 0.5 * DZITRP(I,1)
4218    c  20 CONTINUE
4219          DO L =1,NLEV
4220           DO I =1,irun
4221            IF ( INTCHG(I,L).EQ.1 ) THEN
4222             DDD(I,L) = (         CP  *(DTHV(I,L+1)*DPK(I,L)
4223         .                            + DTHV(I,L)*DPK(I,L+1)) )
4224         .            - ( (2.*RICR)  * ( DU(I,L+1)* DU(I,L)
4225         .                            +   DV(I,L+1)* DV(I,L)) )
4226             AAA(I,L) = STBFCN(I,L) + STBFCN(I,L+1)
4227             BBB(I,L) = STBFCN(I,L) - STBFCN(I,L+1)
4228             CCC(I,L) =            1. / BBB(I,L)
4229             DZITRP(I,L) = AAA(I,L) * CCC(I,L)
4230             AAA(I,L) = AAA(I,L) - DDD(I,L)
4231             DDD(I,L) =        (   DDD(I,L) *    DDD(I,L) )
4232         .            - 4. * (STBFCN(I,L+1) * STBFCN(I,L) )
4233             DDD(I,L) = DDD(I,L)*CCC(I,L)*CCC(I,L)
4234             DDD(I,L) =    SQRT( DDD(I,L) )
4235            ENDIF
4236  C  C
4237        IF( INTCHG(I,1).EQ.1 .AND. AAA(I,1).NE.0. ) THEN          IF( INTCHG(I,L).EQ.1 .AND. AAA(I,L).NE.0. ) THEN
4238        DZITRP(I,1) = ( BBB(I,1)*(1.-DDD(I,1)) ) / AAA(I,1)           DZITRP(I,L) = ( BBB(I,L)*(1.-DDD(I,L)) ) / AAA(I,L)
4239        ENDIF          ENDIF
4240  C  C
4241        DZITRP(I,1) = 0.5 * DZITRP(I,1)          DZITRP(I,L) = 0.5 * DZITRP(I,L)
4242     20 CONTINUE         ENDDO
4243          ENDDO
4244  C  C
4245        RETURN        RETURN
4246        END        END
# Line 3785  C*************************************** Line 4274  C***************************************
4274    
4275  C Argument List Declarations  C Argument List Declarations
4276        integer nlev,nlay,irun        integer nlev,nlay,irun
4277        real RI(irun,NLEV),STRT(irun,NLEV),DW2(irun,NLEV)  c     _RL RI(irun,NLEV),STRT(irun,NLEV),DW2(irun,NLEV)
4278        real XL(irun,NLEV),ZKM(irun,NLEV),ZKH(irun,NLEV)  c     _RL XL(irun,NLEV),ZKM(irun,NLEV),ZKH(irun,NLEV)
4279        real QE(irun,NLEV),QQE(irun,NLEV)  c     _RL QE(irun,NLEV),QQE(irun,NLEV)
4280        INTEGER INTSTB(irun,nlev)  c     INTEGER INTSTB(irun,nlev)
4281        real EE(irun,nlay-1),RF(irun,nlay-1)  c     _RL EE(irun,nlay-1),RF(irun,nlay-1)
4282          _RL RI(irun*NLEV,1),STRT(irun*NLEV,1),DW2(irun*NLEV,1)
4283          _RL XL(irun*NLEV,1),ZKM(irun*NLEV,1), ZKH(irun*NLEV,1)
4284          _RL QE(irun*NLEV,1),QQE(irun*NLEV,1)
4285          INTEGER INTSTB(irun*NLEV,1)
4286          _RL EE(irun*(nlay-1),1),RF(irun*(nlay-1),1)
4287    
4288  C Local Variables  C Local Variables
4289        real b1,b2,d3,rf1,rf2,d3b2,d2,e5,d4,d1,d1half,d2half        _RL b1,b2,d3,rf1,rf2,d3b2,d2,e5,d4,d1,d1half,d2half
4290        real rfc,ricr,ch,cm        _RL rfc,ricr,ch,cm
4291        PARAMETER ( B1      =   16.6    )          PARAMETER ( B1      =   16.6    )
4292        PARAMETER ( B2      =   10.1    )        PARAMETER ( B2      =   10.1    )
4293        PARAMETER ( D3      = 0.29397643 )        PARAMETER ( D3      = 0.29397643 )
4294        PARAMETER ( RF1     = 0.2340678 )        PARAMETER ( RF1     = 0.2340678 )
4295        PARAMETER ( RF2     = 0.2231172 )        PARAMETER ( RF2     = 0.2231172 )
4296        PARAMETER ( D3B2    = D3 / RF1  )        PARAMETER ( D3B2    = D3 / RF1  )
4297        PARAMETER ( D2      = RF1       )        PARAMETER ( D2      = RF1       )
4298        PARAMETER ( E5      = 49.66     )          PARAMETER ( E5      = 49.66     )
4299        PARAMETER ( D4      = 2.6532122E-2 )        PARAMETER ( D4      = 2.6532122E-2 )
4300        PARAMETER ( D1      = D4 * E5   )        PARAMETER ( D1      = D4 * E5   )
4301        PARAMETER ( D1HALF = 0.5 * D1 )        PARAMETER ( D1HALF = 0.5 * D1 )
# Line 3888  C*************************************** Line 4382  C***************************************
4382    
4383  C Argument list Declarations  C Argument list Declarations
4384        integer nlev,nlay,irun        integer nlev,nlay,irun
4385        real Q(irun,NLEV),XL(irun,NLEV),STRT(irun,NLEV)  c     _RL Q(irun,NLEV),XL(irun,NLEV),STRT(irun,NLEV)
4386        real DW2(irun,NLEV)  c     _RL DW2(irun,NLEV)
4387        INTEGER INTSTB(irun,nlay), INTQ(irun,nlay)  c     INTEGER INTSTB(irun,nlay), INTQ(irun,nlay)
4388        real ZKM(irun,NLEV),ZKH(irun,NLEV),P3(irun,NLEV)  c     _RL ZKM(irun,NLEV),ZKH(irun,NLEV),P3(irun,NLEV)
4389          _RL Q(irun*NLEV,1),XL(irun*NLEV,1),STRT(irun*NLEV,1)
4390          _RL DW2(irun*NLEV,1)
4391          INTEGER INTSTB(irun*nlay,1), INTQ(irun*nlay,1)
4392          _RL ZKM(irun*NLEV,1),ZKH(irun*NLEV,1),P3(irun*NLEV,1)
4393    
4394  C Local Variables  C Local Variables
4395        real a1,a2,a4,c1,a5,a3,b1,b2,b3,ff2,ff3,ff4        _RL a1,a2,a4,c1,a5,a3,b1,b2,b3,ff2,ff3,ff4
4396        PARAMETER ( A1      =   0.92    )        PARAMETER ( A1      =   0.92    )
4397        PARAMETER ( A2      =   0.74    )        PARAMETER ( A2      =   0.74    )
4398        PARAMETER ( A4      = 6. * A1 * A1)        PARAMETER ( A4      = 6. * A1 * A1)
4399        PARAMETER ( C1      =   0.08    )        PARAMETER ( C1      =   0.08    )
4400        PARAMETER ( A5      = 3.*C1*(-1.) )        PARAMETER ( A5      = 3.*C1*(-1.) )
4401        PARAMETER ( A3      = A4 * A5*(-1.) )        PARAMETER ( A3      = A4 * A5*(-1.) )
4402        PARAMETER ( B1      =   16.6    )          PARAMETER ( B1      =   16.6    )
4403        PARAMETER ( B2      =   10.1    )        PARAMETER ( B2      =   10.1    )
4404        PARAMETER ( B3      = 1. / B1  )          PARAMETER ( B3      = 1. / B1  )
4405        PARAMETER ( FF2     = 9. * A1 * A2 )        PARAMETER ( FF2     = 9. * A1 * A2 )
4406        PARAMETER ( FF3     = (3.*A2*B2) - (9.*A2*A2 )  )        PARAMETER ( FF3     = (3.*A2*B2) - (9.*A2*A2 )  )
4407        PARAMETER ( FF4     = (3.*A2*B2) + (12.*A1*A2 )  )        PARAMETER ( FF4     = (3.*A2*B2) + (12.*A1*A2 )  )
4408    
4409        real F2(irun,nlay-1),F3(irun,nlay-1)        _RL F2(irun*(nlay-1),1),F3(irun*(nlay-1),1)
4410        real F4(irun,nlay-1),XQ(irun,nlay-1)        _RL F4(irun*(nlay-1),1),XQ(irun*(nlay-1),1)
4411    
4412        integer istnlv        integer istnlv
4413        integer i        integer i
# Line 3987  C*************************************** Line 4485  C***************************************
4485    
4486  C Argument List Declarations  C Argument List Declarations
4487        integer nlev,itype,irun        integer nlev,itype,irun
4488        real XX1(irun,NLEV+1),XX2(irun,NLEV+1)        _RL XX1(irun,NLEV+1),XX2(irun,NLEV+1)
4489        real RHOKDZ(irun,NLEV),FLXFAC(irun,NLEV+1)        _RL RHOKDZ(irun,NLEV),FLXFAC(irun,NLEV+1)
4490        real DXX1G(irun),DXX2G(irun)        _RL DXX1G(irun),DXX2G(irun)
4491        real epsl        _RL epsl
4492    C
4493          _RL AA(irun,nlev), BB(irun,nlev), CC(irun,nlev+1)
4494    c     integer istnlv,istnm1,nlevp1,istnlx
4495          integer NLEVP1,NLEVX
4496          integer I,L
4497  C  C
4498        real AA(irun,nlev), BB(irun,nlev), CC(irun,nlev+1)  c     ISTNLV = irun * NLEV
4499        integer istnlv,istnm1,nlevp1,istnlx  c     ISTNM1 = ISTNLV - irun
       integer i  
 C  
       ISTNLV = irun * NLEV  
       ISTNM1 = ISTNLV - irun  
4500        NLEVP1 = NLEV + 1        NLEVP1 = NLEV + 1
4501        ISTNLX = ISTNM1  c     ISTNLX = ISTNM1
4502        IF(ITYPE.EQ.2) ISTNLX = ISTNLV  c     IF(ITYPE.EQ.2) ISTNLX = ISTNLV
4503          NLEVX =  NLEV - 1
4504          IF(ITYPE.EQ.2) NLEVX =  NLEV
4505  C  C
4506  C   DEFINE MATRIX  C   DEFINE MATRIX
4507  C  C
4508        DO 10 I=1,irun        DO 10 I=1,irun
4509        CC(I,1) = 0.        CC(I,1) = 0.
4510     10 CONTINUE     10 CONTINUE
4511        DO 20 I=1,ISTNLX  c     DO 20 I=1,ISTNLX
4512        CC(I,2) = RHOKDZ(I,1) * FLXFAC(I,2)  c     CC(I,2) = RHOKDZ(I,1) * FLXFAC(I,2)
4513     20 CONTINUE  c  20 CONTINUE
4514        DO 30 I=1,ISTNLV        DO L =1,NLEVX
4515        BB(I,1) =   RHOKDZ(I,1) * FLXFAC(I,1)         DO I =1,irun
4516        AA(I,1) = 1. + CC(I,1) + BB(I,1)          CC(I,L+1) = RHOKDZ(I,L) * FLXFAC(I,L+1)
4517     30 CONTINUE         ENDDO
4518          ENDDO
4519    c     DO 30 I=1,ISTNLV
4520    c     BB(I,1) =   RHOKDZ(I,1) * FLXFAC(I,1)
4521    c     AA(I,1) = 1. + CC(I,1) + BB(I,1)
4522    c  30 CONTINUE
4523          DO L =1,NLEV
4524           DO I =1,irun
4525            BB(I,L) =   RHOKDZ(I,L) * FLXFAC(I,L)
4526            AA(I,L) = 1. + CC(I,L) + BB(I,L)
4527           ENDDO
4528          ENDDO
4529  C  C
4530  C     ADD IMPLICIT BACKWARD FORCING FOR QQ  C     ADD IMPLICIT BACKWARD FORCING FOR QQ
4531        IF(ITYPE.EQ.1) THEN        IF(ITYPE.EQ.1) THEN
4532        DO 40 I=1,ISTNLV  c     DO 40 I=1,ISTNLV
4533        AA(I,1) = AA(I,1) - XX2(I,1)  c     AA(I,1) = AA(I,1) - XX2(I,1)
4534     40 CONTINUE  c  40 CONTINUE
4535           DO L =1,NLEV
4536            DO I =1,irun
4537             AA(I,L) = AA(I,L) - XX2(I,L)
4538            ENDDO
4539           ENDDO
4540        ENDIF        ENDIF
4541  C  C
4542  C     SOLVE MATRIX EQUATION FOR XX1  C     SOLVE MATRIX EQUATION FOR XX1
# Line 4045  C Line 4562  C
4562  C  C
4563  C     ELIMINATE UNDERFLOW  C     ELIMINATE UNDERFLOW
4564        IF(ITYPE.EQ.1) THEN        IF(ITYPE.EQ.1) THEN
4565        DO 70 I=1,ISTNLV  c     DO 70 I=1,ISTNLV
4566        IF( XX1(I,1).LT.EPSL ) XX1(I,1) = 0.  c     IF( XX1(I,1).LT.EPSL ) XX1(I,1) = 0.
4567     70 CONTINUE  c  70 CONTINUE
4568           DO L =1,NLEV
4569            DO I =1,irun
4570             IF ( XX1(I,L).LT.EPSL ) XX1(I,L) = 0.
4571            ENDDO
4572           ENDDO
4573        ENDIF        ENDIF
4574  C  C
4575        RETURN        RETURN
# Line 4056  C Line 4578  C
4578        implicit none        implicit none
4579    
4580        integer k,irun        integer k,irun
4581        real A(irun,K),B(irun,K),C(irun,K),Y(irun,K+1)        _RL A(irun,K),B(irun,K),C(irun,K),Y(irun,K+1)
4582        real F(irun,K)        _RL F(irun,K)
4583    
4584        integer i,L,Lm1        integer i,L,Lm1
4585  C  C
# Line 4087  C Line 4609  C
4609        implicit none        implicit none
4610    
4611        integer k,irun        integer k,irun
4612        real A(irun,K),B(irun,K),Y(irun,K+1)        _RL A(irun,K),B(irun,K),Y(irun,K+1)
4613    
4614        integer i,L        integer i,L
4615  C  C
# Line 4104  C Line 4626  C
4626        implicit none        implicit none
4627    
4628        integer k,irun        integer k,irun
4629        real A(irun,K),B(irun,K),C(irun,K),F(irun,K)        _RL A(irun,K),B(irun,K),C(irun,K),F(irun,K)
4630        real Y(irun,K+1)        _RL Y(irun,K+1)
4631    
4632        integer i,L        integer i,L
4633  C  C
# Line 4172  C*************************************** Line 4694  C***************************************
4694    
4695  C Argument List Declarations  C Argument List Declarations
4696        integer nn,irun,itype        integer nn,irun,itype
4697        real VRIB1(IRUN),VRIB2(IRUN)        _RL VRIB1(IRUN),VRIB2(IRUN)
4698        real VWS1(IRUN),VWS2(IRUN),VZ1(IRUN),VUSTAR(IRUN)        _RL VWS1(IRUN),VWS2(IRUN),VZ1(IRUN),VUSTAR(IRUN)
4699        integer IWATER(IRUN)        integer IWATER(IRUN)
4700        real VAPSIM(IRUN),VAPSIHG(IRUN)        _RL VAPSIM(IRUN),VAPSIHG(IRUN)
4701        real VPSIH(IRUN),VPSIG(IRUN),VX(IRUN)        _RL VPSIH(IRUN),VPSIG(IRUN),VX(IRUN)
4702        real VX0(IRUN),VY(IRUN),VY0(IRUN)        _RL VX0(IRUN),VY(IRUN),VY0(IRUN)
4703        LOGICAL LWATER        LOGICAL LWATER
4704        real VDZETA(IRUN),VDZ0(IRUN),VDPSIM(IRUN)        _RL VDZETA(IRUN),VDZ0(IRUN),VDPSIM(IRUN)
4705        real VDPSIH(IRUN)        _RL VDPSIH(IRUN)
4706        integer INTRIB(IRUN)        integer INTRIB(IRUN)
4707        real VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)        _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
4708        real VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)        _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
4709        real VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)        _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
4710        real VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)        _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
4711        real VXNUM(irun),VDZETA1(irun),VDZETA2(irun)        _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
4712        real VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)        _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
4713        real VDPSIMC(irun),VDPSIHC(irun)        _RL VDPSIMC(irun),VDPSIHC(irun)
4714    
4715  C Local Variables  C Local Variables
4716        real xx0max,prfac,xpfac,difsqt,ustz0s,h0byz0,usth0s        _RL xx0max,prfac,xpfac,difsqt,ustz0s,h0byz0,usth0s
4717        PARAMETER ( XX0MAX  =   1.49821 )        PARAMETER ( XX0MAX  =   1.49821 )
4718        PARAMETER ( PRFAC  = 0.595864   )        PARAMETER ( PRFAC  = 0.595864   )
4719        PARAMETER ( XPFAC  = .55        )          PARAMETER ( XPFAC  = .55        )
4720        PARAMETER ( DIFSQT  = 3.872983E-3)        PARAMETER ( DIFSQT  = 3.872983E-3)
4721        PARAMETER ( USTZ0S =   0.2030325E-5)        PARAMETER ( USTZ0S =   0.2030325E-5)
4722        PARAMETER ( H0BYZ0 =    30.0    )        PARAMETER ( H0BYZ0 =    30.0    )
4723        PARAMETER ( USTH0S =  H0BYZ0*USTZ0S )        PARAMETER ( USTH0S =  H0BYZ0*USTZ0S )
4724    
4725        integer VINT1(irun),VINT2(irun)        integer VINT1(irun),VINT2(irun)
4726        real getcon,vk,bmdl,b2uhs        _RL getcon,vk,bmdl,b2uhs
4727        integer i        integer i
4728  C  C
4729        vk = getcon('VON KARMAN')        vk = getcon('VON KARMAN')
# Line 4460  C        COMPUTE ROUGHNESS LENGTH FOR OC Line 4982  C        COMPUTE ROUGHNESS LENGTH FOR OC
4982  C          BASED ON FUNCTIONS OF LARGE AND POND  C          BASED ON FUNCTIONS OF LARGE AND POND
4983  C          AND OF KONDO --- DESIGNED FOR K = .4  C          AND OF KONDO --- DESIGNED FOR K = .4
4984  C *********************************************************************  C *********************************************************************
4985        implicit none        implicit none
4986    
4987  C Argument List Delcarations  C Argument List Delcarations
4988        integer irun        integer irun
4989        real VZSEA(IRUN),VUSTAR(IRUN),VDZSEA(IRUN)        _RL VZSEA(IRUN),VUSTAR(IRUN),VDZSEA(IRUN)
4990        integer IWATER(IRUN)        integer IWATER(IRUN)
4991        LOGICAL LDZSEA        LOGICAL LDZSEA
4992    
4993  C Local Variables  C Local Variables
4994        real USTMX1,USTMX2,USTMX3        _RL USTMX1,USTMX2,USTMX3
4995        PARAMETER ( USTMX1 =   1.14973  )        PARAMETER ( USTMX1 =   1.14973  )
4996        PARAMETER ( USTMX2 =   0.381844 )        PARAMETER ( USTMX2 =   0.381844 )
4997        PARAMETER ( USTMX3 =   0.0632456)        PARAMETER ( USTMX3 =   0.0632456)
4998    
4999        real AA(IRUN,5),TEMP(IRUN)        _RL AA(IRUN,5),TEMP(IRUN)
5000        integer INT2(IRUN),INT3(IRUN),INT4(IRUN)        integer INT2(IRUN),INT3(IRUN),INT4(IRUN)
5001        integer i,k        integer i,k
5002    
5003        real AA1(5),AA2(5),AA3(5),AA4(5)        _RL AA1(5),AA2(5),AA3(5),AA4(5)
5004        DATA AA1/.2030325E-5,0.0,0.0,0.0,0.0/        DATA AA1/.2030325E-5,0.0,0.0,0.0,0.0/
5005        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,
5006       1         0.395649E-04/       1         0.395649E-04/
# Line 4563  C Line 5085  C
5085        END        END
5086    
5087        subroutine seaice ( nocean, timstp, hice,        subroutine seaice ( nocean, timstp, hice,
5088       .                     eturb,  dedtc,         .                     eturb,  dedtc,
5089       .                    hsturb, dhsdtc,       .                    hsturb, dhsdtc,
5090       .                      qice,  dqice,       .                      qice,  dqice,
5091       .                     swnet,  lwnet, dst4,       .                     swnet,  lwnet, dst4,
5092       .                       pke,  seaic, tc, qa )       .                       pke,  seaic, tc, qa )
5093        implicit none        implicit none
5094        integer  nocean        integer  nocean
5095        real timstp        _RL timstp
5096        real eturb(nocean),dedtc(nocean),hsturb(nocean),dhsdtc(nocean)        _RL eturb(nocean),dedtc(nocean),hsturb(nocean),dhsdtc(nocean)
5097        real swnet(nocean),lwnet(nocean),  dst4(nocean)        _RL swnet(nocean),lwnet(nocean),  dst4(nocean)
5098        real  qice(nocean),dqice(nocean)        _RL  qice(nocean),dqice(nocean)
5099        real   pke(nocean),   tc(nocean),    qa(nocean)        _RL   pke(nocean),   tc(nocean),    qa(nocean)
5100        real seaic(nocean)        _RL seaic(nocean)
5101    
5102  C  rho*C = 1.93e6 J/(m**3 K) ; Peixoto & Oort  C  rho*C = 1.93e6 J/(m**3 K) ; Peixoto & Oort
5103        real rhoC        _RL rhoC
5104        parameter (rhoC = 1.93e6)        parameter (rhoC = 1.93e6)
5105    
5106        real faceps,getcon,latent,codt,deltg,hice        _RL faceps,getcon,latent,codt,deltg,hice
5107        integer i        integer i
5108    
5109        faceps = getcon('EPSFAC')        faceps = getcon('EPSFAC')
5110        latent = getcon('HEATI') * getcon('CALTOJ')        latent = getcon('HEATI') * getcon('CALTOJ')
5111  C Note hice is in centimeters  C Note hice is in centimeters
5112        codt   = rhoC * (hice/100) / timstp          codt   = rhoC * (hice/100) / timstp
5113    
5114  c Update TC and QA  C Update TC and QA
5115  c ----------------  C ----------------
5116        do i =1,nocean        do i =1,nocean
5117         if( seaic(i).gt.0.0 ) then         if( seaic(i).gt.0.0 ) then
5118          deltg = ( swnet(i)-lwnet(i)-latent*eturb(i)-hsturb(i)+qice(i) )          deltg = ( swnet(i)-lwnet(i)-latent*eturb(i)-hsturb(i)+qice(i) )
5119       .          / ( codt+dst4(i)+latent*dedtc(i)+dhsdtc(i)-dqice(i) )       .          / ( codt+dst4(i)+latent*dedtc(i)+dhsdtc(i)-dqice(i) )
5120          qa(i) = qa(i) + (faceps*qa(i)/(tc(i)*tc(i)))*deltg          qa(i) = qa(i) + (faceps*qa(i)/(tc(i)*tc(i)))*deltg
5121          tc(i) = tc(i) + deltg          tc(i) = tc(i) + deltg
5122         endif         endif
5123        enddo        enddo
5124    
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.47
  ViewVC Help
Powered by ViewVC 1.1.22