File: C:\NOAA\NEMS_11731\src\atmos\phys\gbphys.f

1     !===================================================================== !
2     !  description:                                                         !
3     !                                                                       !
4     !     gbphys is the driver subroutine to invoke GFS AM physics          !
5     !     (except radiation but radiative heating is applied here)          !
6     !     at physics time steps                                             !
7     !                                                                       !
8     !  usage:                                                               !
9     !                                                                       !
10     !    call gbphys                                                        !
11     !       inputs:                                                         !
12     !         ( im,ix,levs,lsoil,lsm,ntrac,ncld,ntoz,ntcw,                  !
13     !           nmtvr,nrcm,ko3,lonf,latg,jcap,num_p3d,num_p2d,              !
14     !           kdt,lat,me,pl_coeff,nlons,ncw,flgmin,crtrh,cdmbgwd,         !
15     !           ccwf,dlqf,ctei_rm,clstp,dtp,dtf,fhour,solhr,                !
16     !           slag,sdec,cdec,sinlat,coslat,pgr,ugrs,vgrs,                 !
17     !           tgrs,qgrs,vvel,prsi,prsl,prslk,prsik,phii,phil,             !
18     !           rann,prdout,poz,dpshc,hprime,xlon,xlat,                     !
19     !           slope,shdmin,shdmax,snoalb,tg3,slmsk,vfrac,                 !
20     !           vtype,stype,uustar,oro,coszen,sfcdsw,sfcnsw,                !
21     !           sfcdlw,tsflw,sfcemis,sfalb,swh,hlw,ras,pre_rad,             !
22     !           ldiag3d,lggfs3d,lgocart,lssav,lssav_cc,                     !
23     !           xkzm_m,xkzm_h,xkzm_s,psautco,prautco,evpco,                 !
24     !           flipv,old_monin,cnvgwd,shal_cnv,sashal,newsas,cal_pre,      !
25     !           mom4ice,mstrat,trans_trac,nst_fcst,moist_adj,fscav,         !
26     !           thermodyn_id, sfcpress_id, gen_coord_hybrid,                !
27     !       input/outputs:                                                  !
28     !           hice,fice,tisfc,tsea,tprcp,cv,cvb,cvt,                      !
29     !           srflag,snwdph,weasd,sncovr,zorl,canopy,                     !
30     !           ffmm,ffhh,f10m,srunoff,evbsa,evcwa,snohfa,                  !
31     !           transa,sbsnoa,snowca,soilm,tmpmin,tmpmax,                   !
32     !           dusfc,dvsfc,dtsfc,dqsfc,totprcp,gflux,                      !
33     !           dlwsfc,ulwsfc,suntim,runoff,ep,cldwrk,                      !
34     !           dugwd,dvgwd,psmean,cnvprcp,spfhmin,spfhmax,rain,rainc,      !
35     !           dt3dt,dq3dt,du3dt,dv3dt,dqdt_vacv,acvb,acvt,                !
36     !           slc,smc,stc,upd_mf,dwn_mf,det_mf,dkh,rnp,phy_f3d,phy_f2d,   !
37     !           dlwsfc_cc,ulwsfc_cc,dtsfc_cc,swsfc_cc,                      !
38     !           dusfc_cc, dvsfc_cc, dqsfc_cc,precr_cc,                      !
39     !           xt,xs,xu,xv,xz,zm,xtts,xzts,d_conv,ifd,dt_cool,Qrain,      !
40     !       outputs:                                                        !
41     !           gt0,gq0,gu0,gv0,t2m,q2m,u10m,v10m,                          !
42     !           zlvl,psurf,hpbl,pwat,t1,q1,u1,v1,                           !
43     !           chh,cmm,dlwsfci,ulwsfci,dswsfci,uswsfci,                    !
44     !           dtsfci,dqsfci,gfluxi,epi,smcwlt2,smcref2,wet1,              !
45     !hchuang code change
46     !           gsoil,gtmp2m,gustar,gpblh,gu10m,gv10m,gzorl,goro,           !
47     !           xmu_cc,dlw_cc,dsw_cc,snw_cc,lprec_cc,                       !
48     !           tref, z_c, c_0, c_d, w_0, w_d, rqtk )                       !
49     !                                                                       !
50     !  subprograms called:                                                  !
51     !                                                                       !
52     !     get_prs,  dcyc2t2_pre_rad (testing),    dcyc2t3,  sfc_diff,       !
53     !     sfc_ocean,sfc_drv,  sfc_land, sfc_sice, sfc_diag, moninp1,        !
54     !     moninp,   moninq1,  moninq,   gwdps,    ozphys,   get_phi,        !
55     !     sascnv,   sascnvn,  rascnv,   gwdc,     shalcvt3, shalcv,         !
56     !     shalcnv,  cnvc90,   lrgscl,   gsmdrive, gscond,   precpd,         !
57     !     progt2.                                                           !
58     !                                                                       !
59     !                                                                       !
60     !  program history log:                                                 !
61     !           19xx  - ncep mrf/gfs                                        !
62     !           2002  - s. moorthi  modify and restructure and add Ferrier  !
63     !                               microphysics as an option               !
64     !           200x  - h. juang    modify (what?)                                  !
65     !      nov  2004  - x. wu       modify sea-ice model                    !
66     !      may  2005  - s. moorthi  modify and restructure                  !
67     !           2005  - s. lu       modify to include noah lsm              !
68     !      oct  2006  - h. wei      modify lsm options to include both      !
69     !                               noah and osu lsms.                      !
70     !           2006  - s. moorthi  added a. johansson's convective gravity !
71     !                               wave parameterization code              !
72     !           2007  - s. moorthi  added j. han's modified pbl/sas options !
73     !      dec  2007  - xu li       modified the operational version for    !
74     !                               nst model                               !
75     !           2008  - s. moorthi  applied xu li's nst model to new gfs    !
76     !      mar  2008  - y.-t. hou   added sunshine duration var (suntim) as !
77     !                     an input/output argument.                         !
78     !           2008  - jun wang    added spfhmax/spfhmin as input/output.  !
79     !      apr  2008  - y.-t. hou   added lw sfc emissivity var (sfcemis),  !
80     !                     define the lw sfc dn/up fluxes in two forms: atmos!
81     !                     and ground. also changed sw sfc net flux direction!
82     !                     (positive) from ground -> atmos to the direction  !
83     !                     of atmos -> ground. recode the program and add    !
84     !                     program documentation block.
85     !           2008/ - s. moorthi and y.t. hou upgraded the code to more   !
86     !           2009    modern form and changed all the inputs to MKS units.!
87     !      feb  2009  - s. moorthi  upgraded to add Hochun's gocart changes !
88     !      jul  2009  - s. moorthi  added rqtk for sela's semi-lagrangian   !
89     !      aug  2009  - s. moorthi  added j. han and h. pan updated shallow !
90     !                               convection package                      !
91     !      sep  2009  - s. moorthi  updated for the mcica (rrtm3) radiation !
92     !      dec  2010  - sarah lu    lgocart added to input arg;             !
93     !                               compute dqdt_v if inline gocart is on   !
94     !                                                                       !
95     !                                                                       !
96     !  ====================  defination of variables  ====================  !
97     !                                                                       !
98     !  inputs:                                                       size   !
99     !     ix, im   - integer, horiz dimention and num of used pts      1    !
100     !     levs     - integer, vertical layer dimension                 1    !
101     !     lsoil    - integer, number of soil layers                    1    !
102     !     lsm      - integer, flag for land surface model to use       1    !
103     !                =0  for osu lsm; =1  for noah lsm                      !
104     !     ntrac    - integer, number of tracers                        1    !
105     !     ncld     - integer, number of cloud species                  1    !
106     !     ntoz     - integer, ozone location in the tracer array       1    !
107     !     ntcw     - integer, cloud condensate location in the tracer  1    !
108     !                         array                                    1    !
109     !     nmtvr    - integer, number of topographic variables such as  1    !
110     !                         variance etc used in the GWD parameterization !
111     !     nrcm     - integer, second dimension for the random number   1    !
112     !                         array rann                                    !
113     !     ko3      - integer, number of layers for ozone data          1    !
114     !     lonf,latg- integer, number of lon/lat points                 1    !
115     !     jcap     - integer, number of spectral wave trancation       1    !
116     !                         used only by sascnv shalcnv                   !
117     !     num_p3d  - integer, number of 3D arrays needed for           1    !
118     !                          microphysics                                 !
119     !     num_p2d  - integer, number of 2D arrays needed for           1    !
120     !                         microphysics                                  !
121     !     kdt       -integer, number of the current time step          1    !
122     !     lat       -integer, latitude index - used for debug prints   1    !
123     !     me        -integer, pe number - used for debug prints        1    !
124     !     pl_coeff - integer, number coefficients in ozone forcing     1    !
125     !     nlons    - integer, number of total grid points in a latitude     !
126     !                         circle through a point                   im   !
127     !     ncw      - integer, range of droplet number concentrations for    !
128     !                         Ferrier microphysics                     2    !
129     !     flgmin   - real, range of  minimum large ice fraction for         !
130     !                         Ferrier microphys                        2    !
131     !     crtrh    - real, critical relative humidity at the surface, PBL   !
132     !                      top and at the top of the atmosphere        3    !
133     !     cdmbgwd  - real, multiplication factors for cdmb and gwd     2    !
134     !     ccwf     - real, multiplication factor for critical cloud         !
135     !                      workfunction for RAS                        2    !
136     !     dlqf     - real, factor for cloud condensate detrainment from     !
137     !                      cloud edges (RAS)                           2    !
138     !     ctei_rm  - real, critical cloud top entrainment instability  2    !
139     !                      criteria (used if mstrat=.true.)                 !
140     !     clstp    - real, index used by cnvc90 (for convective clouds)1    !
141     !                      legacy stuff - does not affect forecast          !
142     !     dtp      - real, physics time step in seconds                1    !
143     !     dtf      - real, dynamics time step in seconds               1    !
144     !     fhour    - real, forecast hour                               1    !
145     !     solhr    - real, fcst hour at the end of prev time step      1    !
146     !     slag     - real, equation of time ( radian )                 1    !
147     !     sdec,cdec- real, sin and cos of the solar declination angle  1    !
148     !     sinlat   - real, sin of latitude                             im   !
149     !     coslat   - real, cos of latitude                             im   !
150     !     pgr      - real, surface pressure (Pa)                       im   !
151     !     ugrs,vgrs- real, u/v component of layer wind              ix,levs !
152     !     tgrs     - real, layer mean temperature ( k )             ix,levs !
153     !     qgrs     - real, layer mean tracer concentration     ix,levs,ntrac!
154     !     vvel     - real, layer mean vertical velocity (Pa/s)      ix,levs !
155     !     prsi     - real, pressure at layer interfaces             ix,levs+1
156     !     prsl     - real, mean layer presure                       ix,levs !
157     !     prsik    - real, Exner function at layer interface        ix,levs+1
158     !     prslk    - real, Exner function at layer                  ix,levs !
159     !     phii     - real, interface geopotential height            ix,levs+1
160     !     phil     - real, layer geopotential height                ix,levs !
161     !     rann     - real, random number array (0-1)                ix,nrcm !
162     !     prdout   - real, ozone forcing data                       ix,ko3,pl_coeff!
163     !     poz      - real, ozone forcing data level pressure (ln(Pa))  ko3  !
164     !     dpshc    - real, maximum pressure depth for shallow convection im   !
165     !     hprime   - real, orographic std dev                       ix,nmtvr!
166     !     xlon,xlat- real, longitude and latitude ( radian )           im   !
167     !     slope    - real, sfc slope type for lsm                      im   !
168     !     shdmin   - real, min fractional coverage of green veg        im   !
169     !     shdmax   - real, max fractnl cover of green veg (not used)   im   !
170     !     snoalb   - real, max snow albedo over land (for deep snow)   im   !
171     !     tg3      - real, deep soil temperature                       im   !
172     !     slmsk    - real, sea/land/ice mask (=0/1/2)                  im   !
173     !     vfrac    - real, vegetation fraction                         im   !
174     !     vtype    - real, vegetation type                             im   !
175     !     stype    - real, soil type                                   im   !
176     !     uustar   - real, boundary layer parameter                    im   !
177     !     oro      - real, orography                                   im   !
178     !     coszen   - real, avg cosz over daytime sw radiation interval im   !
179     !     sfcdsw   - real, total sky sfc downward sw flux ( w/m**2 )   im   !
180     !     sfcnsw   - real, total sky sfc netsw flx into ground(w/m**2) im   !
181     !     sfcdlw   - real, total sky sfc downward lw flux ( w/m**2 )   im   !
182     !     tsflw    - real, sfc air (layer 1) temp over lw interval (k) im   !
183     !     sfcemis  - real, sfc lw emissivity ( fraction )              im   !
184     !     sfalb    - real, mean sfc diffused sw albedo                 im   !
185     !     swh      - real, total sky sw heating rates ( k/s )       ix,levs !
186     !     hlw      - real, total sky lw heating rates ( k/s )       ix,levs !
187     !     ras      - logical, flag for ras convection scheme           1    !
188     !     pre_rad  - logical, flag for testing purpose                 1    !
189     !     ldiag3d  - logical, flag for 3d diagnostic fields            1    !
190     !     lggfs3d  - logical, flag for 3d diagnostic fields for gocart 1    !
191     !     lgocart  - logical, flag for 3d diagnostic fields for gocart 1    !
192     !     lssav    - logical, flag controls data store and output      1    !
193     !     lssav_cc - logical, flag for save data for ocean coupling    1    !
194     !     flipv    - logical, flag for vertical direction flip (ras)   1    !
195     !     xkzm_m   - real, background vertical diffusion for momentum  1    !
196     !     xkzm_h   - real, background vertical diffusion for heat, q   1    !
197     !     xkzm_s   - real, sigma threshold for background mom. diffusn 1    !
198     !     psautco  - real, auto conversion coeff from ice to snow      2    !
199     !     prautco  - real, auto conversion coeff from cloud to rain    2    !
200     !     evpco    - real, coeff for evaporation of largescale rain    1    !
201     !     old_monin- logical, flag for diff monin schemes              1    !
202     !     cnvgwd   - logical, flag for conv gravity wave drag          1    !
203     !     shal_cnv - logical, flag for calling shallow convection      1    !
204     !     sashal   - logical, flag for new shallow conv scheme         1    !
205     !     newsas   - logical, flag for new sas conv scheme             1    !
206     !     cal_pre  - logical, flag controls precip type algorithm      1    !
207     !     mom4ice  - logical, flag controls mom4 sea-ice               1    !
208     !     mstrat   - logical, flag for moorthi approach for stratus    1    !
209     !     trans_trac-logical, flag for convective transport of tracers 1    !
210     !     nst_fcst  -integer, flag 0 for no nst, 1 for uncoupled nst        !
211     !                          and 2 for coupled NST                   1    !
212     !     moist_adj- logical, flag for moist convective adjustment     1    !
213     !     fscav    - real, tracer convective scavenging coefficient ntrac-ncld-1!
214     !     thermodyn_id - integer, valid for GFS only for get_prs/phi   1    !
215     !     sfcpress_id  - integer, valid for GFS only for get_prs/phi   1    !
216     !     gen_coord_hybrid - logical for Henry's gen coord             1    !
217     !                                                                       !
218     !  input/outputs:                                                       !
219     !     hice     - real, sea-ice thickness                           im   !
220     !     fice     - real, sea-ice concentration                       im   !
221     !     tisfc    - real, sea-ice temperature                         im   !
222     !     tsea     - real, ground surface temperature ( k )            im   !
223     !     tprcp    - real, total precipitation                         im   !
224     !     the following three variables do not affect the forecast          !
225     !     cv,       -real, convective clouds amountt                   im   !
226     !     cvb       -real, convective clouds base pressure (kPa)       im   !
227     !     cvt       -real, convective clouds top  pressure (kPa)       im   !
228     !     srflag   - real, snow/rain flag for precipitation            im   !
229     !     snwdph   - real, actual snow depth (mm) over land/sea ice    im   !
230     !     weasd    - real, water equiv of accumulated  snow depth (kg/m**2)
231     !                      over land and sea ice                       im   !
232     !     sncovr   - real, snow cover over land                        im   !
233     !     zorl     - real, surface roughness                           im   !
234     !     canopy   - real, canopy water                                im   !
235     !     ffmm     - real, fm parameter from PBL scheme                im   !
236     !     ffhh     - real, fh parameter from PBL scheme                im   !
237     !     f10m     - real, fm at 10m                                   im   !
238     !     srunoff  - real, surface water runoff (from lsm)             im   !
239     !     evbsa    - real, noah lsm diagnostics                        im   !
240     !     evcwa    - real, noah lsm diagnostics                        im   !
241     !     snohfa   - real, noah lsm diagnostics                        im   !
242     !     transa   - real, noah lsm diagnostics                        im   !
243     !     sbsnoa   - real, noah lsm diagnostics                        im   !
244     !     snowca   - real, noah lsm diagnostics                        im   !
245     !     soilm    - real, soil moisture                               im   !
246     !     tmpmin   - real, min temperature at 2m height (k)            im   !
247     !     tmpmax   - real, max temperature at 2m height (k)            im   !
248     !     dusfc    - real, u component of surface stress               im   !
249     !     dvsfc    - real, v component of surface stress               im   !
250     !     dtsfc    - real, sensible heat flux (w/m2)                   im   !
251     !     dqsfc    - real, latent heat flux (w/m2)                     im   !
252     !     totprcp  - real, accumulated total precipitation (kg/m2)     im   !
253     !     gflux    - real, groud conductive heat flux                  im   !
254     !     dlwsfc   - real, time accumulated sfc dn lw flux ( w/m**2 )  im   !
255     !     ulwsfc   - real, time accumulated sfc up lw flux ( w/m**2 )  im   !
256     !     suntim   - real, sunshine duration time (s)                  im   !
257     !     runoff   - real, total water runoff                          im   !
258     !     ep       - real, potential evaporation                       im   !
259     !     cldwrk   - real, cloud workfunction (valid only with sas)    im   !
260     !     dugwd    - real, vertically integrated u change by OGWDi     im   !
261     !     dvgwd    - real, vertically integrated v change by OGWD      im   !
262     !     psmean   - real, surface pressure (kPa)                      im   !
263     !     cnvprcp  - real, accumulated convective precipitation (kg/m2)im   !
264     !     spfhmin  - real, minimum specific humidity                   im   !
265     !     spfhmax  - real, maximum specific humidity                   im   !
266     !     rain     - real, total rain at this time step                im   !
267     !     rainc    - real, convective rain at this time step           im   !
268     !     dt3dt    - real, temperature change due to physics           ix,levs,6 !
269     !     dq3dt    - real, moisture change due to physics              ix,levs,5+pl_coeff!
270     !     du3dt    - real, u momentum change due to physics            ix,levs,4 !
271     !     dv3dt    - real, v momentum change due to physics            ix,levs,4 !
272     !     dqdt_v   - real, total moisture tendency (kg/kg/s)           ix,levs   !
273     !     acv      - real,  array containing accumulated convective clouds im   !
274     !     acvb,acvt- real,  arrays used by cnvc90                      im   !
275     !     slc      - real, liquid soil moisture                     ix,lsoil!
276     !     smc      - real, total soil moisture                      ix,lsoil!
277     !     stc      - real, soil temperature                         ix,lsoil!
278     !     upd_mf   - real, convective updraft mass flux             ix,levs !
279     !     dwn_mf   - real, convective downdraft mass flux           ix,levs !
280     !     det_mf   - real, convective detrainment mass flux         ix,levs !
281     !     dkh      - real, vertical diffusion coefficient (gocart)  ix,levs !
282     !     rnp      - real, n cloud precipitation rate     (gocart)  ix,levs !
283     !     phy_f3d  - real, 3d arrays saved for restart              ix,levs,num_p3d!
284     !     phy_f2d  - real, 2d arrays save for restart               ix,num_p2d!
285     !     dlwsfc_cc- real, sfc dnwd lw flux (w/m**2) for ocn coupling  im   !
286     !     ulwsfc_cc- real, sfc upwd lw flux (w/m**2) for ocn coupling  im   !
287     !     swsfc_cc - real, sfc net sw  flux (w/m**2) for ocn coupling  im   !
288     !     dusfc_cc - real, sfc u-wind                for ocn coupling  im   !
289     !     dvsfc_cc - real, sfc v-wind                for ocn coupling  im   !
290     !     dtsfc_cc - real, sfc temperature  (k)      for ocn coupling  im   !
291     !     dqsfc_cc - real, sfc pressure              for ocn coupling  im   !
292     !     precr_cc - real, total precipitation       for ocn coupling  im   !
293     !
294     !     xt       - real, heat content in DTL                         im   !
295     !     xs       - real, salinity  content in DTL                    im   !
296     !     xu       - real, u-current content in DTL                    im   !
297     !     xv       - real, v-current content in DTL                    im   !
298     !     xz       - real, DTL thickness                               im   !
299     !     zm       - real, MXL thickness                               im   !
300     !     xtts     - real, d(xt)/d(ts)                                 im   !
301     !     xzts     - real, d(xz)/d(ts)                                 im   !
302     !     d_conv   - real, thickness of Free Convection Layer (FCL)    im   !
303     !     ifd      - real, index to start DTM run or not               im   !
304     !     dt_cool  - real, Sub-layer cooling amount                    im   !
305     !     Qrain    - real, sensible heat flux due to rainfall (watts)  im   !
306     !                                                                       !
307     !  outputs:                                                             !
308     !     gt0      - real, updated temperature                        ix,levs !
309     !     gq0      - real, updated tracers                            ix,levs,ntrac!
310     !     gu0      - real, updated zonal wind                         ix,levs !
311     !     gv0      - real, update meridional wind                     ix,levs !
312     !     t2m,q2m  - real, 2 meter temperature and humidity            im   !
313     !     u10m,v10m- real, 10 meater u/v wind speed                    im   !
314     !     zlvl     - real, layer 1 height (m)                          im   !
315     !     psurf    - real, surface pressure (Pa)                       im   !
316     !     hpbl     - real, pbl height (m)                              im   !
317     !     pwat     - real, precipitable water                          im   !
318     !     t1       - real, layer 1 temperature (K)                     im   !
319     !     q1       - real, layer 1 specific humidity (kg/kg)           im   !
320     !     u1       - real, layer 1 zonal wind (m/s)                    im   !
321     !     v1       - real, layer 1 merdional wind (m/s)                im   !
322     !     chh      - real, thermal exchange coefficient                im   !
323     !     cmm      - real, momentum exchange coefficient               im   !
324     !     dlwsfci  - real, instantaneous sfc dnwd lw flux ( w/m**2 )   im   !
325     !     ulwsfci  - real, instantaneous sfc upwd lw flux ( w/m**2 )   im   !
326     !     dswsfci  - real, instantaneous sfc dnwd sw flux ( w/m**2 )   im   !
327     !     uswsfci  - real, instantaneous sfc upwd sw flux ( w/m**2 )   im   !
328     !     dtsfci   - real, instantaneous sfc sensible heat flux        im   !
329     !     dqsfci   - real, instantaneous sfc latent heat flux          im   !
330     !     gfluxi   - real, instantaneous sfc ground heat flux          im   !
331     !     epi      - real, instantaneous sfc potential evaporation     im   !
332     !     smcwlt2  - real, wilting point (volumetric)                  im   !
333     !     smcref2  - real, soil moisture threshold (volumetric)        im   !
334     !     wet1     - real, normalized soil wetness (for GOCART)        im   !
335     !
336     !     gsoil    - real                                              im   !
337     !     gtmp2m   - real                                              im   !
338     !     gustar   - real                                              im   !
339     !     gpblh    - real                                              im   !
340     !     gu10m    - real                                              im   !
341     !     gv10m    - real                                              im   !
342     !     gzorl    - real                                              im   !
343     !     goro     - real                                              im   !
344     !   
345     !     xmu_cc   - real, cosine of zenith angle at time step         im   !
346     !     dlw_cc   - real, sfc dnwd lw flux at time step for ocn cpl   im   !
347     !     dsw_cc   - real, sfc dnwd sw flux at time step for ocn cpl   im   !
348     !     snw_cc   - real, lower atms snow fall rate for ocn cpl       im   !
349     !     lprec_cc - real, lower atms rain fall rate for ocn cpl       im   !
350     
351     !     tref     - real, Reference Temperature                       im   !
352     !     z_c      - real, Sub-layer cooling thickness                 im   !
353     !     c_0      - real, coefficient1 to calculate d(Tz)/d(Ts)       im   !
354     !     c_d      - real, coefficient2 to calculate d(Tz)/d(Ts)       im   !
355     !     w_0      - real, coefficient3 to calculate d(Tz)/d(Ts)       im   !
356     !     w_d      - real, coefficient4 to calculate d(Tz)/d(Ts)       im   !
357     !     rqtk     - real, mass change due to moisture variation       im   !
358     !                                                                       !
359     !                                                                       !
360     !  ====================    end of description    =====================  !
361     
362           subroutine gbphys                                                 &
363     !  ---  inputs:
364          &    ( im,ix,levs,lsoil,lsm,ntrac,ncld,ntoz,ntcw,                  &
365          &      nmtvr,nrcm,ko3,lonf,latg,jcap,num_p3d,num_p2d,              &
366          &      kdt,lat,me,pl_coeff,nlons,ncw,flgmin,crtrh,cdmbgwd,         &
367          &      ccwf,dlqf,ctei_rm,clstp,dtp,dtf,fhour,solhr,                &
368          &      slag,sdec,cdec,sinlat,coslat,pgr,ugrs,vgrs,                 &
369          &      tgrs,qgrs,vvel,prsi,prsl,prslk,prsik,phii,phil,             &
370          &      rann,prdout,poz,dpshc,hprime,xlon,xlat,                     &
371          &      slope,shdmin,shdmax,snoalb,tg3,slmsk,vfrac,                 &
372          &      vtype,stype,uustar,oro,coszen,sfcdsw,sfcnsw,                &
373          &      sfcdlw,tsflw,sfcemis,sfalb,swh,hlw,ras,pre_rad,             &
374     !    &      ldiag3d,lggfs3d,lssav,                                      &
375     !    &      ldiag3d,lggfs3d,lssav,lssav_cc,                             &
376          &      ldiag3d,lggfs3d,lgocart,lssav,lssav_cc,                     &
377          &      xkzm_m,xkzm_h,xkzm_s,psautco,prautco,evpco,                 &
378          &      flipv,old_monin,cnvgwd,shal_cnv,sashal,newsas,cal_pre,      &
379          &      mom4ice,mstrat,trans_trac,nst_fcst,moist_adj,fscav,         &
380          &      thermodyn_id, sfcpress_id, gen_coord_hybrid,                &
381     !  ---  input/outputs:
382          &      hice,fice,tisfc,tsea,tprcp,cv,cvb,cvt,                      &
383          &      srflag,snwdph,weasd,sncovr,zorl,canopy,                     &
384          &      ffmm,ffhh,f10m,srunoff,evbsa,evcwa,snohfa,                  &
385          &      transa,sbsnoa,snowca,soilm,tmpmin,tmpmax,                   &
386          &      dusfc,dvsfc,dtsfc,dqsfc,totprcp,gflux,                      &
387          &      dlwsfc,ulwsfc,suntim,runoff,ep,cldwrk,                      &
388          &      dugwd,dvgwd,psmean,cnvprcp,spfhmin,spfhmax,rain,rainc,      &
389          &      dt3dt,dq3dt,du3dt,dv3dt,dqdt_v,acv,acvb,acvt,               &
390          &      slc,smc,stc,upd_mf,dwn_mf,det_mf,dkh,rnp,phy_f3d,phy_f2d,   &
391          &      dlwsfc_cc,ulwsfc_cc,dtsfc_cc,swsfc_cc,                      &
392          &      dusfc_cc, dvsfc_cc, dqsfc_cc,precr_cc,                      &
393          &      xt,xs,xu,xv,xz,zm,xtts,xzts,d_conv,ifd,dt_cool,Qrain,       &
394     !  ---  outputs:
395          &      gt0,gq0,gu0,gv0,t2m,q2m,u10m,v10m,                          &
396          &      zlvl,psurf,hpbl,pwat,t1,q1,u1,v1,                           &
397          &      chh,cmm,dlwsfci,ulwsfci,dswsfci,uswsfci,                    &
398          &      dtsfci,dqsfci,gfluxi,epi,smcwlt2,smcref2,wet1,              &
399          &      gsoil,gtmp2m,gustar,gpblh,gu10m,gv10m,gzorl,goro,           &
400          &      xmu_cc,dlw_cc,dsw_cc,snw_cc,lprec_cc,                       &
401          &      tref, z_c, c_0, c_d, w_0, w_d,                              &
402          &      rqtk                                                        &
403          &      )
404     
405     !
406           use machine ,   only : kind_phys
407           use physcons,   only : con_cp, con_fvirt, con_g, con_rd, con_rv,  &
408          &                       con_hvap, con_hfus, con_rerth, con_pi
409     !
410           implicit none
411     !
412     
413     !  ---  constant parameters:
414           real(kind=kind_phys), parameter :: hocp    = con_hvap/con_cp
415           real(kind=kind_phys), parameter :: fhourpr = 0.0
416           real(kind=kind_phys), parameter :: rlapse  = 0.65e-2
417           real(kind=kind_phys), parameter :: rhc_max = 0.9999
418           real(kind=kind_phys), parameter :: qmin    = 1.0e-10
419           real(kind=kind_phys), parameter :: p850    = 85000.0
420           real(kind=kind_phys), parameter :: epsq    = 1.e-20
421           real(kind=kind_phys), parameter :: hsub    = con_hvap+con_hfus
422     !     real(kind=kind_phys), parameter :: cb2mb   = 10.0
423           real(kind=kind_phys), parameter :: pa2mb   = 0.01
424           real(kind=kind_phys), parameter :: czmin   = 0.0001      ! cos(89.994)
425     
426           real(kind=kind_phys), parameter ::                                &
427     !    & dxmax=ln(1.0/(5000.0*2500.0)),  dxmin=ln(1.0/(192.0*94.0)
428     !    & dxmax=-16.118095651,dxmin=-9.800790154,dxinv=1.0/(dxmax-dxmin)
429          & dxmaxs=-16.118095651,dxmins=-9.800790154,                        &
430          & dxinvs=1.0/(dxmaxs-dxmins),
431     !    & dxmax=ln(1.0/(2000.0*1000.0)),  dxmin=ln(1.0/(192.0*94.0)
432          & dxmaxr=-14.5086577385,dxminr=-9.800790154,                       &
433          & dxinvr=1.0/(dxmaxr-dxminr)
434     !    & dxmax=ln(1.0/(2500.0*1250.0)),  dxmin=ln(1.0/(192.0*94.0)
435     !    & dxmax=-14.95494484115,dxmin=-9.800790154, dxinv=1.0/(dxmax-dxmin)
436            real(kind=kind_phys) dxmax, dxmin, dxinv
437     
438     !  ---  inputs:
439     !  note: lgocart is the logical for in-line gocart;
440     !        lggfs3d is the logical for off-line gocart 
441           integer, intent(in) :: ix,   im,   levs, lsoil,   lsm,     ntrac, &
442          &                       ncld, ntoz, ntcw, nmtvr,   nrcm,    ko3,   &
443          &                       lonf, latg, jcap, num_p3d, num_p2d, kdt,   &
444          &                       me,   pl_coeff, lat,                       &
445          &                       thermodyn_id, sfcpress_id
446     
447     
448           integer, intent(in) :: nlons(im), ncw(2), nst_fcst
449     
450           logical, intent(in) :: ras,        pre_rad,   ldiag3d, flipv,     &
451          &                       old_monin,  cnvgwd,    sashal,  newsas,    &
452          &                       lssav,      lssav_cc,  mom4ice, mstrat,    &
453          &                       trans_trac, moist_adj, lggfs3d, cal_pre,   &
454          &                       shal_cnv, gen_coord_hybrid, lgocart
455     
456           real(kind=kind_phys), dimension(im),            intent(in) ::     &
457          &      sinlat, coslat, pgr,    dpshc,  xlon,   xlat,               &
458          &      slope,  shdmin, shdmax, snoalb, tg3,    slmsk,  vfrac,      &
459          &      vtype,  stype,  uustar, oro,    coszen, sfcnsw, sfcdsw,     &
460          &      sfcdlw, tsflw,  sfalb,  sfcemis
461     
462           real(kind=kind_phys), dimension(ix,levs),       intent(in) ::     &
463          &      ugrs, vgrs, tgrs, vvel, prsl, prslk, phil, swh, hlw
464     
465           real(kind=kind_phys), intent(in) ::  qgrs(ix,levs,ntrac)
466     
467           real(kind=kind_phys), dimension(ix,levs+1),     intent(in) ::     &
468          &      prsi, prsik, phii
469     
470           real(kind=kind_phys), intent(in) ::  hprime(ix,nmtvr),            &
471          &      prdout(ix,ko3,pl_coeff),       rann(ix,nrcm), poz(ko3)
472     
473           real(kind=kind_phys), intent(in) ::  dtp,     dtf, fhour, solhr,  &
474          &      slag,    sdec,     cdec,       ctei_rm(2), clstp,           &
475          &      ccwf(2), crtrh(3), flgmin(2),  dlqf(2), cdmbgwd(2),         &
476          &      xkzm_m,  xkzm_h, xkzm_s, psautco(2),   prautco(2), evpco
477     
478           real(kind=kind_phys), intent(in) ::  fscav(ntrac-ncld-1)
479     
480     
481     !  ---  input/output:
482           real(kind=kind_phys), dimension(im),            intent(inout) ::  &
483          &      hice,   fice,    tisfc,  tsea,   tprcp,  cv,     cvb,  cvt, &
484          &      srflag, snwdph,  weasd, sncovr, zorl,   canopy, ffmm, ffhh, &
485          &      f10m,   srunoff, evbsa,  evcwa,  snohfa, transa, sbsnoa,    &
486          &      snowca, soilm,   tmpmin, tmpmax, dusfc,  dvsfc,  dtsfc,     &
487          &      dqsfc,  totprcp, gflux,  dlwsfc, ulwsfc, suntim, runoff, ep,&
488          &      cldwrk, dugwd,   dvgwd,  psmean, cnvprcp,spfhmin, spfhmax,  &
489          &      rain, rainc, acv,    acvb,    acvt,                         &
490     ! om coupling
491          &      dlwsfc_cc, ulwsfc_cc, dtsfc_cc, swsfc_cc, dusfc_cc,         &
492          &      dvsfc_cc,  dqsfc_cc,  precr_cc,                             &
493     ! nst
494          &      xt(im),     xs(im),   xu(im),     xv(im),  xz(im),   zm(im),&
495          &      xtts(im),   xzts(im), d_conv(im), ifd(im), dt_cool(im),     &
496          &      Qrain(im)
497     !
498           real(kind=kind_phys), dimension(ix,lsoil),      intent(inout) ::  &
499          &      smc, stc, slc
500     
501           real(kind=kind_phys), dimension(ix,levs),       intent(inout) ::  &
502          &      upd_mf, dwn_mf, det_mf, dkh, rnp
503     
504           real(kind=kind_phys),                           intent(inout) ::  &
505          &      phy_f3d(ix,levs,num_p3d), phy_f2d(ix,num_p2d),              &
506          &      dt3dt(ix,levs,6), du3dt(ix,levs,4), dv3dt(ix,levs,4),       &
507          &      dq3dt(ix,levs,5+pl_coeff)
508     !    &      dq3dt(ix,levs,5+pl_coeff), dqdt_v(ix,levs)
509     
510     !  ---  output:
511           real(kind=kind_phys), dimension(im),            intent(out) ::    &
512          &      t2m,     q2m,     u10m,    v10m,    zlvl,   psurf, hpbl,    &
513          &      pwat,    t1,      q1,      u1,      v1,     chh,   cmm,     &
514          &      dlwsfci, ulwsfci, dswsfci, uswsfci, dtsfci, dqsfci,         &
515          &      gfluxi,  epi,     smcwlt2, smcref2, wet1,                   &
516     !hchuang code change 11/12/2007 [+2L]
517          &      gsoil(im), gtmp2m(im), gustar(im), gpblh(im), gu10m(im),    &
518          &      gv10m(im), gzorl(im), goro(im),                             &
519          &      xmu_cc,  dlw_cc,  dsw_cc,  snw_cc,  lprec_cc,               &
520          &      tref,    z_c,     c_0,     c_d,     w_0,   w_d, rqtk
521     
522           real(kind=kind_phys), dimension(ix,levs),       intent(out) ::    &
523          &      gt0, gu0, gv0, dqdt_v
524     !    &      gt0, gu0, gv0
525     
526           real(kind=kind_phys), dimension(ix,levs,ntrac), intent(out) ::    &
527          &      gq0
528     
529     !  ---  local:
530           real(kind=kind_phys), dimension(im)          :: ccwfac, garea,    &
531          &      dlength, xncw,   cumabs, qmax,   cice,    zice,   tice,     &
532          &      gflx,                    rainl,  rain1,   raincs, evapc,    &
533     !    &      gflx,    rain,   rainc,  rainl,  rain1,   raincs, evapc,    &
534          &      snowmt,  cd,     cdq,    qss,    dusfcg,  dvsfcg, dusfc1,   &
535          &      dvsfc1,  dtsfc1, dqsfc1, rb,     rhscnpy, drain,  cld1d,    &
536          &      evap,    hflx,   stress, t850,   ep1d,    gamt,   gamq,     &
537          &      sigmaf,                  oc,     theta,   gamma,  sigma,    &
538          &      elvmax,  wind,   work1,  work2,  runof,   xmu,    oro_land, &
539          &      fm10,    fh2,    tsurf,  tx1,    tx2,     ctei_r, flgmin_l, &
540          &      evbs,    evcw,   trans,  sbsno,  snowc,   adjsfcdsw,        &
541          &      adjsfcnsw, adjsfcdlw, adjsfculw, asfcdlw, asfculw, gsfcdlw, &
542          &      gsfculw, xcosz,  tseal,  snohf,  dlqfac,  work3,            &
543          &      domr,    domzr,  domip,  doms,   psautco_l, prautco_l,      &
544          &      ctei_rml
545     
546     !    &      dswsfc, radsl,                                              &
547     !    &      dlwsf1,  ulwsf1, xcosz,  tseal,  snohf,   dlqfac,           &
548     !    &      domr,    domzr,  domip,  doms
549     
550     !     real(kind=kind_phys), dimension(ix,levs)     :: ud_mf, dd_mf,     &
551     !    &      dt_mf, del
552           real(kind=kind_phys), dimension(ix,levs)     :: del
553           real(kind=kind_phys), dimension(im,levs-1)   :: dkt
554     
555           real(kind=kind_phys), dimension(im,levs)     :: rhc, sr, dtdt,    &
556          &      dudt, dvdt, gwdcu, gwdcv, diagn1, diagn2, cuhr, cumchr,     &
557          &      qr_col, fc_ice, rainp, ud_mf, dd_mf, dt_mf
558     
559           real(kind=kind_phys), dimension(im,lsoil)    :: smsoil, stsoil,   &
560          &      ai, bi, cci, rhsmc, zsoil, slsoil
561     
562           real(kind=kind_phys) :: rhbbot, rhbtop, rhpbl, frain, f_rain,     &
563          &      f_ice, qi, qw, qr, wc, tem, tem1, tem2,  sume,  sumr, sumq, &
564          &      dqdt(im,levs,ntrac), oa4(im,4), clx(im,4)
565     
566     !           in clw, the first two varaibles are cloud water and ice. 
567     !           from third to ntrac are convective transportable tracers,
568     !           third being the ozone, when ntrac=3 (valid only with ras)
569     
570           real(kind=kind_phys), allocatable :: clw(:,:,:)
571     
572           integer, dimension(im) :: kbot, ktop, kcnv, soiltyp, vegtype,     &
573          &          kpbl, slopetyp, kinver, lmh, levshc
574     
575           integer :: i, nvdiff, kk, ic, k, n, ipr, lv, k1, iter, levshcm,   &
576          &           tracers, trc_shft, tottracer
577     
578           logical, dimension(im) :: flag_iter, flag_guess, invrsn
579     
580           logical :: lprnt
581     
582     !
583     !===> ...  begin here
584     !
585     !     if (ras) then
586     !       dxmax = dxmaxr ; dxmin = dxminr ; dxinv = dxinvr
587     !     else
588     !       dxmax = dxmaxs ; dxmin = dxmins ; dxinv = dxinvs
589     !     endif
590             dxmax = dxmaxs ; dxmin = dxmins ; dxinv = dxinvs
591      
592     !  --- ...  set up check print point (for debugging)
593     !
594     !*************************************************************************
595     !     lprnt = .true.
596           lprnt = .false.
597     !     lprnt = me == 0
598     
599           ipr = 1
600     
601           if (im > 90) then
602             ipr = 91
603            else
604              lprnt = .false.
605            endif
606     !     lprnt = kdt .gt. 0
607     !     do i = 1, im
608     !       work1(1) = xlon(i) * 57.29578
609     !       if (work1(1) >= 180.0) work1(1) = work1(1) - 360.0
610     !       work2(1) = xlat(i) * 57.29578
611     !       print *,' me=',me,' work1=',work1(1),' work2=',work2(1),' i=',i
612     !       lprnt = kdt > 4320
613     !       lprnt = kdt > 0 .and. abs(work1(1)-110.3) < 0.5                 &
614     !    &                  .and. abs(work2(1)-2.0)   < 0.5
615     !       lprnt = kdt >= 14 .and. lat == 43 
616     !       lprnt = kdt >= 0 .and. abs(xlon(i)*57.29578-143.182) < 0.101    &
617     !    &                   .and. abs(xlat(i)*57.29578-6.235)  < 0.101
618     !       lprnt = kdt >= 0 .and. abs(xlon(i)*57.29578-135.0) < 0.201      &
619     !    &                   .and. abs(xlat(i)*57.29578-10.476)  < 0.201
620     !       lprnt = kdt >= 0 .and. abs(xlon(i)*57.29578-110.3) < 0.201      &
621     !    &                   .and. abs(xlat(i)*57.29578-2.0)   < 0.201
622     !       print *,' i=',i,' xlon=',xlon(i)*57.29578,                      &
623     !    &                  ' xlat=',xlat(i)*57.29578,' i=',i,' me=',me
624     !       if (lprnt) then
625     !         ipr = i
626     !         exit
627     !       endif
628     !     enddo
629     
630     !     lprnt = .false.
631     !     if(lprnt) then
632     !       print *,' im=',im,' ix=',ix,' levs=',levs,' lsoil=',lsoil,      &
633     !    &   ' ntrac=',ntrac,' ntoz=',ntoz,' ntcw=',ntcw,' me=',me,         &
634     !    &   ' xlat=',xlat(ipr),' kdt=',kdt,' slmsk=',slmsk(ipr),           &
635     !    &   ' tsea=',tsea(ipr),                                            &
636     !    &   ' tsea=',tsea(ipr),' tref=',tref(ipr),' dt_cool=',dt_cool(ipr),&
637     !    &   ' dt_warm=',2.0*xt(ipr)/xz(ipr),' nrcm=',nrcm,' xlon=',
638     !    &    xlon(ipr), ' tisfc=',tisfc(ipr),                              &
639     !    &   ' dt_warm=',dt_warm(ipr),' nrcm=',nrcm,' xlon=',xlon(ipr),     &
640     !    &   ' sfalb=',sfalb(ipr),' kdt=',kdt
641     !       print *,' pgr=',pgr(ipr),' kdt=',kdt,' ipr=',ipr
642     !       print *,' ipr=',ipr,' phy_f2d=',phy_f2d(ipr,1:num_p2d)
643     !       print *,' ugrs=',ugrs(ipr,:)
644     !       print *,' vgrs=',vgrs(ipr,:)
645     !       print *,' tgrs=',tgrs(ipr,:),' kdt=',kdt,' ipr=',ipr
646     !    &,  ' xlon=',xlon(ipr),' xlat=',xlat(ipr)
647     !       print *,' qgrs=',qgrs(ipr,:,1)
648     !       print *,' ozg=',qgrs(ipr,:,2)
649     !       print *,' clw=',qgrs(ipr,:,3)
650     !    &,  ' xlon=',xlon(ipr),' xlat=',xlat(ipr)
651     !     endif
652     !
653     !*************************************************************************
654     !
655     !     write (0,*) ' ras=',ras,pre_rad,                                  &
656     !    &      ldiag3d,lggfs3d,lssav,lssav_cc,                             &
657     !    &      xkzm_m,xkzm_h,psautco,evpco,                                &
658     !    &      flipv,old_monin,cnvgwd,sashal,newsas,cal_pre,               &
659     !    &      mom4ice,mstrat,trans_trac,nst_fcst,moist_adj
660     
661           do i=1,im
662             lmh(i) = levs
663           enddo
664           nvdiff = ntrac    ! vertical diffusion of all tracers!
665     !
666     !  --- ...  figure out how many extra tracers are there
667     !
668           tottracer = 0            ! no convective transport of tracers
669           if (trans_trac) then
670             if (ntcw > 0) then
671               if (ntoz < ntcw) then
672                 trc_shft = ntcw + ncld - 1
673               else
674                 trc_shft = ntoz
675               endif
676             elseif (ntoz > 0) then
677               trc_shft = ntoz
678             else
679               trc_shft = 1
680             endif
681     
682             tracers   = ntrac - trc_shft
683             tottracer = tracers
684             if (ntoz > 0) tottracer = tottracer + 1  ! ozone is added separately
685           endif
686     
687     !     if (lprnt) print *,' trans_trac=',trans_trac,' tottracer=',       &
688     !    &                   tottracer,' trc_shft=',trc_shft,' kdt=',kdt
689     
690           allocate ( clw(ix,levs,tottracer+2) )
691     
692     !
693           call get_prs(im,ix,levs,ntrac,tgrs,qgrs,                          &
694          &             thermodyn_id, sfcpress_id,                           &
695          &             gen_coord_hybrid,                                    &
696          &             prsi,prsik,prsl,prslk,phii,phil,del)
697     !
698     !     if (lprnt) then
699     !       print *,' prsi=',prsi(ipr,:)
700     !       print *,' prsl=',prsl(ipr,:)
701     !       print *,' del=',del(ipr,:)
702     !     endif
703     !
704           rhbbot = crtrh(1)
705           rhpbl  = crtrh(2)
706           rhbtop = crtrh(3)
707     !
708     
709           do i = 1, im
710             levshc(i) = 0
711           enddo
712     
713           do k = 2, levs
714             do i = 1, im
715               if (prsi(i,1)-prsi(i,k) <= dpshc(i)) levshc(i) = k
716             enddo
717           enddo
718     
719           levshcm = 1
720           do i = 1, im
721             levshcm = max(levshcm, levshc(i))
722           enddo
723     
724     !  --- ...  frain=factor for centered difference scheme correction of rain amount.
725     
726           frain = dtf / dtp
727     
728           do i = 1, im
729             soiltyp(i)  = int( stype(i)+0.5 )
730             sigmaf(i)   = max( vfrac(i),0.01 )
731     !       sigmaf(i)   = max( vfrac(i),0.3 )
732             if (lsm == 0) sigmaf(i)   =  0.5 + vfrac(i) * 0.5
733     
734             vegtype(i)  = int( vtype(i)+0.5 )
735             slopetyp(i) = int( slope(i)+0.5 )    !! clu: slope -> slopetyp
736     
737             if (slmsk(i) == 2.0) then
738               soiltyp(i)  = 9
739               vegtype(i)  = 13
740               slopetyp(i) = 9                    !! clu: qa(slopetyp)
741             endif
742           enddo
743     
744     !  --- ...  transfer soil moisture and temperature from global to local variables
745     
746     !     write(0,*)' in gbphys stc=',stc(1:5,1)
747           do k = 1, lsoil
748             do i = 1, im
749               smsoil(i,k) = smc(i,k)
750               stsoil(i,k) = stc(i,k)
751               slsoil(i,k) = slc(i,k)          !! clu: slc -> slsoil
752             enddo
753           enddo
754     
755     !  --- ...  xw: transfer ice thickness & concentration from global to local variables
756     
757           do i = 1, im
758             zice(i) = hice(i)
759             cice(i) = fice(i)
760             tice(i) = tisfc(i)
761           enddo
762     
763           do k = 1, levs
764             do i = 1, im
765               dudt(i,k) = 0.
766               dvdt(i,k) = 0.
767               dtdt(i,k) = 0.
768             enddo
769           enddo
770     
771           do n = 1, ntrac
772             do k = 1, levs
773               do i = 1, im
774                 dqdt(i,k,n) = 0.
775               enddo
776             enddo
777           enddo
778     
779           do i = 1, im
780             work1(i)    = (log(coslat(i) / (nlons(i)*latg)) - dxmin) * dxinv
781             work1(i)    = max(0.0, min(1.0,work1(i)))
782             work2(i)    = 1.0 - work1(i)
783             psurf(i)    = pgr(i)
784             work3(i)    = prsik(i,1) / prslk(i,1)
785             tem1        = con_rerth * (con_pi+con_pi)*coslat(i)/nlons(i)
786             tem2        = con_rerth * con_pi/latg
787             garea(i)    = tem1 * tem2
788             dlength(i)  = sqrt( tem1*tem1+tem2*tem2 )
789           enddo
790     
791           if (lssav) then
792             do i = 1, im
793               psmean(i) = psmean(i) + pgr(i)*dtf
794             enddo
795           endif
796     
797     !  --- ...  initialize dtdt with heating rate from dcyc2 
798     
799     !     if (lprnt) then
800     !       do ipr=1,im
801     !         print *,' before dcyc2: im=',im,' lsoil=',lsoil,' levs=',levs &
802     !    &,    ' sde=',sdec,' cdec=',cdec,' tsea=',tsea(ipr),' ipr=',ipr    &
803     !    &,    ' lat=',lat,' me=',me,' kdt=',kdt                            &
804     !    &,    ' sfcdlw=',sfcdlw(ipr),' sfcnsw=',sfcnsw(ipr)
805     !         print *,' hlw=',hlw(ipr,:),' me=',me,' lat=',lat,xlon(ipr)
806     !         print *,' swh=',swh(ipr,:),' me=',me,' lat=',lat,xlon(ipr)
807     !       enddo
808     !     endif
809     
810     !  --- ...  adjust mean radiation fluxes and heating rates to fit for
811     !           faster model time steps.
812     !      sw:  using cos of zenith angle as scaling factor
813     !      lw:  using surface air skin temperature as scaling factor
814     
815           if (pre_rad) then
816     
817             call dcyc2t3_pre_rad                                            &
818     !  ---  inputs:
819          &     ( solhr,slag,sdec,cdec,sinlat,coslat,                        &
820          &       xlon,coszen,tsea,tgrs(1,1),tgrs(1,1),                      &
821          &       sfcdsw,sfcnsw,sfcdlw,swh,hlw,                              &
822          &       ix, im, levs,                                              &
823     !  ---  input/output:
824          &       dtdt,                                                      &
825     !  ---  outputs:
826          &       adjsfcdsw,adjsfcnsw,adjsfcdlw,adjsfculw,xmu,xcosz          &
827     !old vars   ( dswsfc,    -radsl,  dlwsf1,   ulwsf1,  xmu,xcosz )
828          &     )
829     
830           else
831     
832             call dcyc2t3                                                    &
833     !  ---  inputs:
834          &     ( solhr,slag,sdec,cdec,sinlat,coslat,                        &
835          &       xlon,coszen,tsea,tgrs(1,1),tsflw,                          &
836          &       sfcdsw,sfcnsw,sfcdlw,swh,hlw,                              &
837          &       ix, im, levs,                                              &
838     !  ---  input/output:
839          &       dtdt,                                                      &
840     !  ---  outputs:
841          &       adjsfcdsw,adjsfcnsw,adjsfcdlw,adjsfculw,xmu,xcosz          &
842     !old vars   ( dswsfc,    -radsl,  dlwsf1,   ulwsf1,  xmu,xcosz )
843          &     )
844     
845           endif
846     
847     !  ---  convert lw fluxes for land/ocean/sea-ice models
848     !  note: adjsfcdlw and adjsfculw are time-step adjusted lw sfc dn/up fluxes.
849     !        those fluxes from dcyc2t3 subprogram are not applied with emissivity
850     !        effect.  one needs to be cautious that, when emis<1, sfc dn/up lw
851     !        fluxes will expressed differently when looking from atmospheric
852     !        (as the output from radiation) to those applied for land/ocean/sea-ice
853     !        model processes. however, the net effects are the same.
854     !
855     !   - flux to/from lnd/oc/ice:  gsfcdlw=emis*adjsfcdlw; gsfculw=emis*adjsfculw
856     !   - flux to/from atmos mdl: asfcdlw=adjsfcdlw; asfculw=emis*adjsfculw+(1-emis)*adjsfcdlw
857             do i = 1, im
858               gsfcdlw(i) = sfcemis(i) * adjsfcdlw(i)                  ! for lnd/ocn/sice
859               gsfculw(i) = sfcemis(i) * adjsfculw(i)                  ! for lnd/ocn/sice
860               asfcdlw(i) = adjsfcdlw(i)                               ! for atmos output
861               asfculw(i) = gsfculw(i) + adjsfcdlw(i) - gsfcdlw(i)     ! for atmos output
862     !org      asfculw(i) = gsfculw(i) + (1.0-sfcemis(i))*adjsfcdlw(i) ! for atmos output
863             enddo
864     
865     !  --- ...  coupling insertion
866     !  note: all radiation fluxes are positive values.
867     !       sfc net sw is defined as sfcnsw = sfcdsw - sfcusw (*** positive downward
868     !       see subr grrad for def)
869     
870           if (lssav_cc) then
871             do i = 1, im
872               xmu_cc(i)    = max( 0.0, min( 1.0, xcosz(i) ))
873               dsw_cc(i)    = adjsfcdsw(i)
874               dlw_cc(i)    = adjsfcdlw(i)
875               dlwsfc_cc(i) = dlwsfc_cc(i) + gsfcdlw(i)
876               ulwsfc_cc(i) = ulwsfc_cc(i) + gsfculw(i)
877               swsfc_cc(i)  = swsfc_cc(i)  - adjsfcnsw(i) ! swsfc is positive upward
878     !         sfcems_cc(i) = sfcemis(i)                  ! for ocn mdl if needed
879             enddo
880           end if
881     
882     !  --- ...  accumulate/save output variables
883          
884           if (lssav) then
885     
886     !  --- ...  sunshine duration time is defined as the length of time (in mdl output
887     !           interval) that solar radiation falling on a plane perpendicular to the
888     !           direction of the sun >= 120 w/m2
889     
890             do i = 1, im
891               if ( xcosz(i) >= czmin ) then   ! zenth angle > 89.994 deg
892                 tem1 = adjsfcdsw(i) / xcosz(i)
893     
894                 if ( tem1 >= 120.0 ) then
895                   suntim(i) = suntim(i) + dtf
896                 endif
897               endif
898             enddo
899     
900     !  --- ...  sfc lw fluxes used by atmospheric model are saved for output
901     
902             do i = 1, im
903               dlwsfc(i) = dlwsfc(i) + asfcdlw(i)*dtf
904               ulwsfc(i) = ulwsfc(i) + asfculw(i)*dtf
905             enddo
906     
907             if (ldiag3d) then
908               do k = 1, levs
909                 do i = 1, im
910                   dt3dt(i,k,1) = dt3dt(i,k,1) + hlw(i,k)*dtf
911                   dt3dt(i,k,2) = dt3dt(i,k,2) + swh(i,k)*dtf*xmu(i)
912                 enddo
913               enddo
914             endif
915     
916           endif    ! end if_lssav_block
917     
918           do i = 1, im
919             kcnv(i)   = 0
920             kinver(i) = levs
921             invrsn(i) = .false.
922             tx1(i)    = 0.0
923             tx2(i)    = 10.0
924             ctei_r(i) = 10.0
925           enddo
926     
927     !    Only used for old shallow convection with mstrat=.true.
928     
929           if (.not. sashal .and. shal_cnv) then
930             do i = 1, im
931               ctei_rml(i) = ctei_rm(1)*work1(i) + ctei_rm(2)*work2(i)
932             enddo
933             do k = 1, levs/2
934               do i = 1, im
935     
936                 if (prsi(i,1)-prsi(i,k+1) < 0.35*prsi(i,1)                  &
937          &          .and. (.not. invrsn(i))) then
938                   tem = (tgrs(i,k+1)-tgrs(i,k)) / (prsl(i,k)-prsl(i,k+1))
939     
940                   if ((tem > 0.00010 .and. tx1(i) < 0.0) .or. 
941          &            (tem-abs(tx1(i)) > 0.0 .and. tx2(i) < 0.0)) then
942                     invrsn(i) = .true.
943     
944                     if (qgrs(i,k,1) > qgrs(i,k+1,1)) then
945                       tem1 = tgrs(i,k+1) + hocp*max(qgrs(i,k+1,1),qmin)
946                       tem2 = tgrs(i,k)   + hocp*max(qgrs(i,k,1),qmin)
947     
948                       tem1 = tem1 / prslk(i,k+1) - tem2 / prslk(i,k)
949     
950     !  --- ...  (cp/l)(deltathetae)/(deltatwater) > ctei_rm -> conditon for CTEI
951                       ctei_r(i) = (1.0/hocp)*tem1/(qgrs(i,k+1,1)-qgrs(i,k,1)&
952          &                      + qgrs(i,k+1,ntcw)-qgrs(i,k,ntcw))
953                     else
954                       ctei_r(i) = 10
955                     endif
956     
957                     if ( ctei_rml(i) > ctei_r(i) ) then
958                       kinver(i) = k
959                     else
960                       kinver(i) = levs
961     !                 kinver(i) = k + 1    ! commented on 12/15/09
962                     endif
963     !               kinver(i) = k
964     !Aug 4          kinver(i) = k + 1
965                   endif
966     
967                   tx2(i) = tx1(i)
968                   tx1(i) = tem
969                 endif
970     
971               enddo
972             enddo
973           endif
974     
975     !  --- ...  check print
976     
977     !     ipr = 1
978     !     if (lprnt) then
979     !       print *,' before progtm: im=',im,' lsoil=',lsoil                &
980     !    &,         ' nvdiff=',nvdiff,' adjsfcnsw=',adjsfcnsw(ipr)          &
981     !    &,         ' asfcdlw=',asfcdlw(ipr),' asfculw=',asfculw(ipr)       &
982     !    &,         ' gsfcdlw=',gsfcdlw(ipr),' gsfculw=',gsfculw(ipr)       &
983     !    &,         ' sfcemis=',sfcemis(ipr),' tsea2=',tsea(ipr)            &
984     !    &,         ' ipr=',ipr,' me=',me,' lat=',lat,' xlon=',xlon(ipr)    &
985     !    &,         ' kdt=',kdt
986     
987     !       print *,' dtdth=',dtdt(ipr,:),' kdt=',kdt
988     !     endif
989     
990     !  --- ...  lu: initialize flag_guess, flag_iter, tsurf
991     
992           do i = 1, im
993             tsurf(i)      = tsea(i)
994             flag_guess(i) = .false.
995             flag_iter(i)  = .true.
996             drain(i)      = 0.0
997             ep1d(i)       = 0.0
998             runof(i)      = 0.0
999             hflx(i)       = 0.0
1000             evap(i)       = 0.0
1001     
1002             evbs(i)       = 0.0
1003             evcw(i)       = 0.0
1004             trans(i)      = 0.0
1005             sbsno(i)      = 0.0
1006             snowc(i)      = 0.0
1007             snohf(i)      = 0.0
1008             zlvl(i)       = phil(i,1) / con_g
1009             smcwlt2(i)    = 0.0
1010             smcref2(i)    = 0.0
1011           enddo
1012     
1013     !  --- ...  lu: iter-loop over (sfc_diff,sfc_drv,sfc_ocean,sfc_sice)
1014     
1015           do iter = 1, 2
1016     
1017     !  --- ...  surface exchange coefficients
1018     !
1019             call sfc_diff(im,pgr,ugrs,vgrs,tgrs,qgrs,zlvl,                  &
1020          &                tsea,zorl,cd,cdq,rb,                              &
1021          &                prsl(1,1),work3,slmsk,                            &
1022          &                stress,ffmm,ffhh,                                 &
1023          +                uustar,wind,phy_f2d(1,num_p2d),fm10,fh2,          &
1024     !   This line below needs to commented for pre-Helin correction to LM - Moorthi
1025          &                sigmaf,vegtype,shdmax,                            &
1026          +                tsurf, flag_iter)
1027     
1028     !       if (lprnt) print *,' cdq=',cdq(ipr),' iter=',iter               &
1029     !    &,   ' wind=',wind(ipr),'phy_f2d=',phy_f2d(ipr,num_p2d),' ugrs='   &
1030     !    &,   ugrs(ipr,1),' vgrs=',vgrs(ipr,1)
1031     
1032     !  --- ...  lu: update flag_guess
1033     
1034             do i = 1, im
1035               if (iter == 1 .and. wind(i) < 2.0) then
1036                 flag_guess(i) = .true.
1037               endif
1038             enddo
1039     
1040     !  --- ...  surface energy balance over ocean
1041     
1042             if ( nst_fcst > 0 ) then
1043     
1044               do i = 1, im
1045                 if ( slmsk(i) == 0.0 ) then
1046                   tseal(i) = tsea(i)  + oro(i) * rlapse
1047                   tsurf(i) = tsurf(i) + oro(i) * rlapse
1048                 endif
1049               enddo
1050               if ( nst_fcst > 1 ) then
1051                 do i = 1, im
1052                   if ( slmsk(i) == 0.0 ) then
1053                     tref(i)  = tseal(i) - (xt(i)+xt(i))/xz(i) + dt_cool(i)
1054                   endif
1055                 enddo
1056               else
1057                 do i = 1, im
1058                   if ( slmsk(i) == 0.0 ) then
1059                     tref(i)  = tseal(i)
1060                   endif
1061                 enddo
1062               endif
1063     
1064     !         if (lprnt) print *,' tseaz1=',tsea(ipr),' tref=',tref(ipr),   &
1065     !    &      ' dt_cool=',dt_cool(ipr),' dt_warm=',2.0*(xt(ipr)/xz(ipr)   &
1066     !    &      ' dt_cool=',dt_cool(ipr),' dt_warm=',dt_warm(ipr)           &
1067     !    &,     ' tgrs=',tgrs(ipr,1),' prsl=',prsl(ipr,1)
1068     !    &,     ' work3=',work3(ipr),' kdt=',kdt
1069     
1070               call sfc_nst                                                  &
1071          &       ( im,lsoil,pgr,ugrs,vgrs,tgrs,qgrs,tref,cd,cdq,            &
1072          &         prsl(1,1),work3,slmsk,xlon,sinlat,stress,                &
1073          &         sfcemis,gsfcdlw,adjsfcnsw,tprcp,dtf,kdt,                 &
1074          &         phy_f2d(1,num_p2d),flag_iter,flag_guess,nst_fcst,        &
1075          &         lprnt,ipr,                                               &
1076     !  --- Input/output
1077          &         tseal,tsurf,xt,xs,xu,xv,xz,zm,xtts,xzts,dt_cool,         &
1078          &         z_c,c_0,c_d,w_0,w_d,d_conv,ifd,Qrain,                    &
1079     
1080     !    &         tseal, ifd, time_old, time_ins, i_sw, i_q,               &
1081     !    &         i_qrain, i_m, i_tau, i_sw_zw, i_q_ts, i_m_ts, dt_cool,   &
1082     !    &         dt_warm, z_c, z_w, c_0, c_d, w_0, w_d,                   &
1083     !  ---  outputs:
1084          &         qss, gflx, cmm, chh, evap, hflx, ep1d)
1085     
1086     !         if (lprnt) print *,' tseaz2=',tseal(ipr),' tref=',tref(ipr),  &
1087     !    &     ' dt_cool=',dt_cool(ipr),' dt_warm=',2.0*xt(ipr)/xz(ipr),    &
1088     !    &     ' kdt=',kdt
1089     !    &     ' dt_cool=',dt_cool(ipr),' dt_warm=',dt_warm(ipr),' kdt=',kdt
1090     
1091               do i = 1, im
1092                 if ( slmsk(i) == 0.0 ) then
1093                   tsurf(i) = tsurf(i) - oro(i) * rlapse
1094                 endif
1095               enddo
1096               if ( nst_fcst > 1 ) then
1097                 do i = 1, im
1098                   if ( slmsk(i) == 0.0 ) then
1099                     tsea(i)  = tref(i) + (xt(i)+xt(i))/xz(i) - dt_cool(i)
1100          &                             - oro(i) * rlapse
1101                   endif
1102                 enddo
1103               endif
1104     
1105     !         if (lprnt) print *,' tseaz2=',tsea(ipr),' tref=',tref(ipr),   &
1106     !    &    ' dt_cool=',dt_cool(ipr),' dt_warm=',dt_warm(ipr),' kdt=',kdt
1107     
1108             else
1109     
1110               call sfc_ocean                                                &
1111     !  ---  inputs:
1112          &     ( im,pgr,ugrs,vgrs,tgrs,qgrs,tsea,cd,cdq,                    &
1113          &       prsl(1,1),work3,slmsk,phy_f2d(1,num_p2d),flag_iter,        &
1114     !  ---  outputs:
1115          &       qss,cmm,chh,gflx,evap,hflx,ep1d                            &
1116          &     )
1117     !
1118             endif
1119     
1120     !       if (lprnt) print *,' sfalb=',sfalb(ipr),' ipr=',ipr             &
1121     !    &,   ' weasd=',weasd(ipr),' snwdph=',snwdph(ipr)                   &
1122     !    &,   ' tprcp=',tprcp(ipr),' kdt=',kdt,' iter=',iter
1123      
1124     !  --- ...  surface energy balance over land
1125     !
1126             if (lsm == 1) then                          ! noah lsm call
1127     
1128               call sfc_drv                                                  &
1129     !  ---  inputs:
1130          &     ( im,lsoil,pgr,ugrs,vgrs,tgrs,qgrs,soiltyp,vegtype,sigmaf,   &
1131          &       sfcemis,gsfcdlw,adjsfcdsw,adjsfcnsw,dtf,tg3,cd,cdq,        &
1132          &       prsl(1,1),work3,zlvl,slmsk,phy_f2d(1,num_p2d),slopetyp,    &
1133          &       shdmin,shdmax,snoalb,sfalb,flag_iter,flag_guess,           &
1134     !  ---  in/outs:
1135          &       weasd,snwdph,tsea,tprcp,srflag,smsoil,stsoil,slsoil,       &
1136          &       canopy,trans,tsurf,                                        &
1137     !  ---  outputs:
1138          &       sncovr,qss,gflx,drain,evap,hflx,ep1d,runof,                &
1139          &       cmm,chh,evbs,evcw,sbsno,snowc,soilm,snohf,                 &
1140          &       smcwlt2,smcref2,wet1                                       &
1141          &     )
1142     
1143             else                                       ! osu lsm call
1144     
1145               call sfc_land                                                 &
1146     !  ---  inputs:
1147          &     ( im,lsoil,pgr,ugrs,vgrs,tgrs,qgrs,smsoil,soiltyp,           &
1148          &       sigmaf,vegtype,sfcemis,gsfcdlw,adjsfcnsw,dtf,              &
1149          &       zorl,tg3,cd,cdq,prsl(1,1),work3,slmsk,                     &
1150          &       phy_f2d(1,num_p2d),flag_iter,flag_guess,                   &
1151     !  ---  input/outputs:
1152          &       weasd,tsea,tprcp,srflag,stsoil,canopy,tsurf,               &
1153     !  ---  outputs:
1154          &       qss,snowmt,gflx,zsoil,rhscnpy,rhsmc,                       &
1155          &       ai,bi,cci,drain,evap,hflx,ep1d,cmm,chh,                    &
1156          &       evbs,evcw,trans,sbsno,snowc,soilm,                         &
1157          &       snohf,smcwlt2,smcref2                                      &
1158          &     )
1159     
1160             endif
1161     
1162     !       if (lprnt) print *,' tseabeficemodel =',tsea(ipr),' me=',me     &
1163     !    &,   ' kdt=',kdt,' stsoil=',stsoil(ipr,:)
1164     
1165     !  --- ...  surface energy balance over seaice
1166     
1167             call sfc_sice                                                   &
1168     !  ---  inputs:
1169          &     ( im,lsoil,pgr,ugrs,vgrs,tgrs,qgrs,dtf,                      &
1170          &       sfcemis,gsfcdlw,adjsfcnsw,adjsfcdsw,srflag,                &
1171          &       cd,cdq,prsl(1,1),work3,slmsk,phy_f2d(1,num_p2d),           &
1172          &       flag_iter,mom4ice,lsm, lprnt,ipr,                          &
1173     !    &       flag_iter,mom4ice,lsm,                                     &
1174     !  ---  input/outputs:
1175          &       zice,cice,tice,weasd,tsea,tprcp,stsoil,ep1d,               &
1176     !  ---  outputs:
1177          &       snwdph,qss,snowmt,gflx,cmm,chh,evap,hflx                   &
1178          &     )
1179     
1180     !  --- ...  lu: update flag_iter and flag_guess
1181     
1182             do i = 1, im
1183               flag_iter(i)  = .false.
1184               flag_guess(i) = .false.
1185     
1186               if(slmsk(i) == 1. .and. iter == 1) then
1187                 if (wind(i) < 2.0) flag_iter(i) = .true.
1188               elseif (slmsk(i) == 0. .and. iter == 1                        &
1189          &                           .and. nst_fcst > 1) then
1190                 if (wind(i) < 2.0) flag_iter(i) = .true.
1191               endif
1192             enddo
1193     
1194           enddo   ! end iter_loop
1195     
1196           do i = 1, im
1197             epi(i)     = ep1d(i)
1198             dlwsfci(i) = gsfcdlw(i)
1199             ulwsfci(i) = gsfculw(i)
1200             uswsfci(i) = adjsfcdsw(i) - adjsfcnsw(i)
1201             dswsfci(i) = adjsfcdsw(i)
1202             gfluxi(i)  = gflx(i)
1203             t1(i)      = tgrs(i,1)
1204             q1(i)      = qgrs(i,1,1)
1205             u1(i)      = ugrs(i,1)
1206             v1(i)      = vgrs(i,1)
1207           enddo
1208     
1209           if (lsm == 0) then                          ! osu lsm call
1210             do i = 1, im
1211              sncovr(i) = 0.0
1212              if (weasd(i) > 0.0) sncovr(i) = 1.0
1213             enddo
1214           endif
1215     
1216     !  --- ...  update near surface fields
1217     
1218           call sfc_diag(im,lsoil,pgr,ugrs,vgrs,tgrs,qgrs,                   &
1219          &              tsea,qss,f10m,u10m,v10m,t2m,q2m,work3,slmsk,        &
1220          &              evap,ffmm,ffhh,fm10,fh2)
1221     
1222           do i = 1, im
1223             phy_f2d(i,num_p2d) = 0.0
1224           enddo
1225     
1226     !     if (lprnt) print *,' tseaim=',tsea(ipr),' me=',me,' kdt=',kdt
1227     
1228           if (lssav) then
1229             do i = 1, im
1230     !     write(0,*)' i=',i,' gflx=',gflx(i),' evbs=',evbs(i),' evcw=',
1231     !    & evcw(i)
1232     !    &,' trans=',trans(i),' sbsno=',sbsno(i),' snowc=',snowc(i),
1233     !    &' shohf=',snohf(i),' i=',i
1234               gflux(i)   = gflux(i)  + gflx(i)*dtf
1235               evbsa(i)   = evbsa(i)  + evbs(i)*dtf
1236               evcwa(i)   = evcwa(i)  + evcw(i)*dtf
1237               transa(i)  = transa(i) + trans(i)*dtf
1238               sbsnoa(i)  = sbsnoa(i) + sbsno(i)*dtf
1239               snowca(i)  = snowca(i) + snowc(i)*dtf
1240               snohfa(i)  = snohfa(i) + snohf(i)*dtf
1241     
1242               tmpmax(i)  = max(tmpmax(i),t2m(i))
1243               tmpmin(i)  = min(tmpmin(i),t2m(i))
1244     
1245               spfhmax(i) = max(spfhmax(i),q2m(i))
1246               spfhmin(i) = min(spfhmin(i),q2m(i))
1247               ep(i)      = ep(i) + ep1d(i) * dtf
1248     !hchuang code change 11/12/2007 [+6]
1249     !hchuang oro is passin variable not derived in this routine
1250               gtmp2m(i)  =  gtmp2m(i) + t2m(i)    * DTF
1251               gu10m(i)   =  gu10m(i)  + u10m(i)   * DTF
1252               gv10m(i)   =  gv10m(i)  + v10m(i)   * DTF
1253               gustar(i)  =  gustar(i) + uustar(i) * DTF
1254               gzorl(i)   =  gzorl(i)  + zorl(i)   * DTF
1255               goro(i)    =  goro(i)   + slmsk(i)  * DTF
1256             enddo
1257           endif
1258     
1259           do i = 1, im
1260     
1261     !  --- ...  compute coefficient of evaporation in evapc
1262     !
1263             if (evapc(i) > 1.0e0) evapc(i) = 1.0e0
1264     
1265     !  --- ...  over snow cover or ice or sea, coef of evap =1.0e0
1266     
1267             if (weasd(i) > 0.0 .or. slmsk(i) /= 1.0) evapc(i) = 1.0e0
1268     
1269           enddo
1270     
1271     !  --- ...  vertical diffusion
1272     
1273     !     if (lprnt) print *,' tsea3=',tsea(ipr),' slmsk=',slmsk(ipr)       &
1274     !    &, ' kdt=',kdt,' evap=',evap(ipr)
1275     !     if (lprnt)  print *,' dtdtb=',dtdt(ipr,:)
1276     
1277     !     do i = 1, im
1278     !       if (slmsk(i) == 0.0) then
1279     !         oro_land(i) = 0.0
1280     !       else
1281     !         oro_land(i) = oro(i)
1282     !       endif
1283     !     enddo
1284     
1285           if (old_monin) then
1286     
1287             if (mstrat) then
1288               call moninp1(ix,im,levs,nvdiff,dvdt,dudt,dtdt,dqdt,           &
1289          &     ugrs,vgrs,tgrs,qgrs,                                         &
1290          &     prsik(1,1),rb,ffmm,ffhh,tsea,qss,hflx,evap,stress,wind,kpbl, &
1291          &     prsi,del,prsl,prslk,phii,phil,dtp,                           &
1292          &     dusfc1,dvsfc1,dtsfc1,dqsfc1,hpbl,gamt,gamq,dkt,              &
1293     !    &     kinver, ctei_r, ctei_rm, xkzm_m, xkzm_h)
1294          &     kinver, xkzm_m, xkzm_h)
1295             else
1296               call moninp(ix,im,levs,nvdiff,dvdt,dudt,dtdt,dqdt,            &
1297          &     ugrs,vgrs,tgrs,qgrs,                                         &
1298          &     prsik(1,1),rb,ffmm,ffhh,tsea,qss,hflx,evap,stress,wind,kpbl, &
1299          &     prsi,del,prsl,prslk,phii,phil,dtp,                           &
1300          &     dusfc1,dvsfc1,dtsfc1,dqsfc1,hpbl,gamt,gamq,dkt,xkzm_m,xkzm_h)
1301             endif
1302     
1303           else
1304     
1305     !       if (mstrat) then
1306     !         call moninq1(ix,im,levs,nvdiff,dvdt,dudt,dtdt,dqdt,           &
1307     !    &     ugrs,vgrs,tgrs,qgrs,swh,hlw,xmu,slmsk,                       &
1308     !    &     prsik(1,1),rb,ffmm,ffhh,tsea,qss,hflx,evap,stress,wind,kpbl, &
1309     !    &     prsi,del,prsl,prslk,phii,phil,dtp,                           &
1310     !    &     dusfc1,dvsfc1,dtsfc1,dqsfc1,hpbl,gamt,gamq,dkt,              &
1311     !    &     kinver)
1312     !       else
1313               call moninq(ix,im,levs,nvdiff,dvdt,dudt,dtdt,dqdt,            &
1314     !    &     ugrs,vgrs,tgrs,qgrs,swh,hlw,xmu,slmsk,                       &
1315          &     ugrs,vgrs,tgrs,qgrs,swh,hlw,xmu,                             &
1316          &     prsik(1,1),rb,ffmm,ffhh,tsea,qss,hflx,evap,stress,wind,kpbl, &
1317          &     prsi,del,prsl,prslk,phii,phil,dtp,                           &
1318          &     dusfc1,dvsfc1,dtsfc1,dqsfc1,hpbl,gamt,gamq,dkt,              &
1319          &     kinver, xkzm_m, xkzm_h, xkzm_s)
1320     !       endif
1321     
1322           endif   ! end if_old_monin
1323     
1324     !     if (lprnt) then
1325     !       print *,' dusfc1=',dusfc1(ipr)
1326     !       print *,' dtsfc1=',dtsfc1(ipr)
1327     !       print *,' dqsfc1=',dqsfc1(ipr)
1328     !       print *,' dtdt=',dtdt(ipr,:)
1329     !       print *,' dudtm=',dudt(ipr,:)
1330     !     endif
1331     
1332     !  --- ...  coupling insertion
1333     
1334           if (lssav_cc) then
1335             do i=1, im
1336               dusfc_cc(i) = dusfc_cc(i) + dusfc1(i)        !*dtf <-na h.(?)
1337               dvsfc_cc(i) = dvsfc_cc(i) + dvsfc1(i)        !*dtf <-na h.(?)
1338               dtsfc_cc(i) = dtsfc_cc(i) + dtsfc1(i)        !*dtf <-na h.(?)
1339               dqsfc_cc(i) = dqsfc_cc(i) + dqsfc1(i)        !*dtf <-na h.(?)
1340             enddo
1341           endif
1342     
1343           if (lssav) then
1344             do i = 1, im
1345               dusfc(i)  = dusfc(i) + dusfc1(i)*dtf
1346               dvsfc(i)  = dvsfc(i) + dvsfc1(i)*dtf
1347               dtsfc(i)  = dtsfc(i) + dtsfc1(i)*dtf
1348               dqsfc(i)  = dqsfc(i) + dqsfc1(i)*dtf
1349               dtsfci(i) = dtsfc1(i)
1350               dqsfci(i) = dqsfc1(i)
1351     !hchuang code change 11/12/2007 [+1L]
1352               gpblh(i)  =  gpblh(i)  + hpbl(i)*dtf
1353             enddo
1354     
1355             if (ldiag3d) then
1356               do k = 1, levs
1357                 do i = 1, im
1358                   tem  = dtdt(i,k) - (hlw(i,k)+swh(i,k)*xmu(i))
1359                   dt3dt(i,k,3) = dt3dt(i,k,3) + tem*dtf
1360     !             dq3dt(i,k,1) = dq3dt(i,k,1) + dqdt(i,k,1) * dtf
1361                   du3dt(i,k,1) = du3dt(i,k,1) + dudt(i,k)   * dtf
1362                   du3dt(i,k,2) = du3dt(i,k,2) - dudt(i,k)   * dtf
1363                   dv3dt(i,k,1) = dv3dt(i,k,1) + dvdt(i,k)   * dtf
1364                   dv3dt(i,k,2) = dv3dt(i,k,2) - dvdt(i,k)   * dtf
1365                 enddo
1366               enddo
1367             endif
1368     ! update dqdt_v to include moisture tendency due to vertical diffusion
1369             if (lgocart) then
1370               do k = 1, levs
1371                 do i = 1, im
1372                   dqdt_v(i,k)  = dqdt(i,k,1) * dtf
1373                 enddo
1374               enddo
1375             endif
1376             if (ldiag3d .or. lggfs3d) then
1377               do k = 1, levs
1378                 do i = 1, im
1379                   tem  = dqdt(i,k,1) * dtf
1380                   dq3dt(i,k,1) = dq3dt(i,k,1) + tem
1381     !             dqdt_v(i,k)  = tem
1382                 enddo
1383               enddo
1384               if (ntoz > 0) then
1385                 do k = 1, levs
1386                   do i = 1, im
1387                     dq3dt(i,k,5) = dq3dt(i,k,5) + dqdt(i,k,ntoz)*dtf
1388                   enddo
1389                 enddo
1390               endif
1391             endif
1392             if (lggfs3d) then
1393               do k=1,levs-1
1394                do i=1,im
1395                   dkh(i,k) = dkh(i,k) + dkt(i,k)*dtf
1396                 enddo
1397               enddo
1398             endif
1399     
1400     
1401           endif   ! end if_lssav
1402     
1403           if (nmtvr == 6) then
1404     
1405             do i = 1, im
1406               oc(i) = hprime(i,2)
1407             enddo
1408     
1409             do k = 1, 4
1410               do i = 1, im
1411                 oa4(i,k) = hprime(i,k+2)
1412                 clx(i,k) = 0.0
1413               enddo
1414             enddo
1415     
1416           elseif (nmtvr == 10) then
1417     
1418             do i = 1, im
1419               oc(i) = hprime(i,2)
1420             enddo
1421     
1422             do k = 1, 4
1423               do i = 1, im
1424                 oa4(i,k) = hprime(i,k+2)
1425                 clx(i,k) = hprime(i,k+6)
1426               enddo
1427             enddo
1428     
1429           elseif (nmtvr == 14) then
1430     
1431             do i = 1, im
1432               oc(i) = hprime(i,2)
1433             enddo
1434     
1435             do k = 1, 4
1436               do i = 1, im
1437                 oa4(i,k) = hprime(i,k+2)
1438                 clx(i,k) = hprime(i,k+6)
1439               enddo
1440             enddo
1441     
1442             do i = 1, im
1443               theta(i)  = hprime(i,11)
1444               gamma(i)  = hprime(i,12)
1445               sigma(i)  = hprime(i,13)
1446               elvmax(i) = hprime(i,14)
1447             enddo
1448     
1449           else
1450     
1451             oc  = 0
1452             oa4 = 0
1453             clx = 0
1454             theta = 0
1455             gamma = 0
1456             sigma = 0
1457             elvmax = 0
1458     
1459           endif   ! end if_nmtvr
1460     
1461           call gwdps(im, ix, im,   levs,  dvdt, dudt,                       &
1462          &           ugrs,   vgrs, tgrs,  qgrs,                             &
1463          &           kpbl,   prsi, del,   prsl, prslk,                      &
1464          &           phii,   phil, dtp,                                     &
1465          &           kdt,    hprime(1,1), oc, oa4, clx,                     &
1466          &           theta,sigma,gamma,elvmax,dusfcg, dvsfcg,               &
1467          &           con_g,con_cp,con_rd,con_rv, lonf, nmtvr, cdmbgwd,      &
1468          &           me, lprnt,ipr)
1469     
1470     !     if (lprnt)  print *,' dudtg=',dudt(ipr,:)
1471     
1472           if (lssav) then
1473             do i = 1, im
1474               dugwd(i) = dugwd(i) + dusfcg(i)*dtf
1475               dvgwd(i) = dvgwd(i) + dvsfcg(i)*dtf
1476             enddo
1477     
1478     !       if (lprnt) print *,' dugwd=',dugwd(ipr),' dusfcg=',dusfcg(ipr)
1479     !       if (lprnt) print *,' dvgwd=',dvgwd(ipr),' dvsfcg=',dvsfcg(ipr)
1480     
1481             if (ldiag3d) then
1482               do k = 1, levs
1483                 do i = 1, im
1484                   du3dt(i,k,2) = du3dt(i,k,2) + dudt(i,k) * dtf
1485                   dv3dt(i,k,2) = dv3dt(i,k,2) + dvdt(i,k) * dtf
1486                 enddo
1487               enddo
1488             endif
1489           endif
1490     
1491           do  k = 1, levs
1492             do i = 1, im
1493               gt0(i,k)   = tgrs(i,k)   + dtdt(i,k)   * dtp
1494               gu0(i,k)   = ugrs(i,k)   + dudt(i,k)   * dtp
1495               gv0(i,k)   = vgrs(i,k)   + dvdt(i,k)   * dtp
1496             enddo
1497           enddo
1498     
1499           do n = 1, ntrac
1500             do k = 1, levs
1501               do i = 1, im
1502                 gq0(i,k,n) = qgrs(i,k,n) + dqdt(i,k,n) * dtp
1503               enddo
1504             enddo
1505           enddo
1506     
1507     !  --- ...  check print
1508     
1509     !     if (me == 0) then
1510     !       sumq = 0.0
1511     !       do k = 1, levs
1512     !         do i = 1, im
1513     !           sumq = sumq + (dqdt(i,k,1)+dqdt(i,k,ntcw)) * del(i,k)
1514     !         enddo
1515     !       enddo
1516     
1517     !       sume = 0.0
1518     !       do i = 1, im
1519     !         sume = sume + dqsfc1(i)
1520     !       enddo
1521     
1522     !       sumq = sumq * 1000.0 / con_g
1523     !       sume = sume / con_hvap
1524     !       print *,' after mon: sumq=',sumq,' sume=',sume, ' kdt=',kdt
1525     !     endif
1526     
1527     !  --- ...  ozone physics
1528     
1529           if (ntoz > 0 .and. ntrac >= ntoz) then
1530     
1531             call ozphys(ix,im,levs,ko3,dtp,gq0(1,1,ntoz),gq0(1,1,ntoz)      &
1532          &,             gt0, poz, prsl, prdout, pl_coeff, del, ldiag3d      &
1533     !hchuang code change [r1L]
1534          &,             lggfs3d,dq3dt(1,1,6), me)
1535     !    &,             dq3dt(1,1,6), me)
1536     
1537           endif
1538     
1539     !  --- ...  to side-step the ozone physics
1540     
1541     !      if (ntrac >= 2) then
1542     !        do k = 1, levs
1543     !          gq0(k,ntoz) = qgrs(k,ntoz)
1544     !        enddo
1545     !      endif
1546     
1547     !     if (lprnt) then
1548     !       print *,' levs=',levs,' jcap=',jcap,' dtp',dtp                  &
1549     !    *,  ' slmsk=',slmsk(ilon,ilat),' kdt=',kdt
1550     !       print *,' rann=',rann,' ncld=',ncld,' iq=',iq,' lat=',lat
1551     !       print *,' pgr=',pgr
1552     !       print *,' del=',del(ipr,:)
1553     !       print *,' prsl=',prsl(ipr,:)
1554     !       print *,' prslk=',prslk(ipr,:)
1555     !       print *,' rann=',rann(ipr,1)
1556     !       print *,' gt0=',gt0(ipr,:)                                      &
1557     !    &,         ' kdt=',kdt,' xlon=',xlon(ipr),' xlat=',xlat(ipr)
1558     !       print *,' dtdt=',dtdt(ipr,:)
1559     !       print *,' gu0=',gu0(ipr,:)
1560     !       print *,' gv0=',gv0(ipr,:)
1561     !       print *,' gq0=',(gq0(ipr,k,1),k=1,levs)
1562     !       print *,' gq1=',(gq0(ipr,k,ntcw),k=1,levs)
1563     !       print *,' vvel=',vvel
1564     !     endif
1565     
1566           if (ldiag3d) then
1567     
1568             do k = 1, levs
1569               do i = 1, im
1570                 dtdt(i,k)   = gt0(i,k)
1571     !           dqdt(i,k,1) = gq0(i,k,1)
1572                 dudt(i,k)   = gu0(i,k)
1573                 dvdt(i,k)   = gv0(i,k)
1574               enddo
1575             enddo
1576     
1577           elseif (cnvgwd) then
1578     
1579             do k = 1, levs
1580               do i = 1, im
1581                 dtdt(i,k)   = gt0(i,k)
1582               enddo
1583             enddo
1584     
1585           endif   ! end if_ldiag3d/cnvgwd
1586           if (ldiag3d .or. lggfs3d) then
1587             do k = 1, levs
1588               do i = 1, im
1589                 dqdt(i,k,1)   = gq0(i,k,1)
1590               enddo
1591             enddo
1592           endif   ! end if_ldiag3d/lggfs3d
1593     
1594           call get_phi(im,ix,levs,ntrac,gt0,gq0,                            &
1595          &             thermodyn_id, sfcpress_id,                           &
1596          &             gen_coord_hybrid,                                    &
1597          &             prsi,prsik,prsl,prslk,phii,phil)
1598     
1599     !     if (lprnt) then
1600     !       print *,' phii2=',phii(ipr,:)
1601     !       print *,' phil2=',phil(ipr,:)
1602     !     endif
1603     
1604           do k = 1, levs
1605             do i = 1, im
1606               clw(i,k,1) = 0.0
1607               clw(i,k,2) = -999.9
1608             enddo
1609           enddo
1610     
1611     !  --- ...  for convective tracer transport (while using ras)
1612     
1613           if (ras) then
1614             if (tottracer > 0) then
1615     
1616               if (ntoz > 0) then
1617                 do k=1,levs
1618                   do i=1,im
1619                     clw(i,k,3) = gq0(i,k,ntoz)
1620                   enddo
1621                 enddo
1622     
1623                 if (tracers > 0) then
1624                   do n = 1, tracers
1625                     do k=1,levs
1626                       do i=1,im
1627                         clw(i,k,3+n) = gq0(i,k,n+trc_shft)
1628                       enddo
1629                     enddo
1630                   enddo
1631                 endif
1632               else
1633                 do n = 1, tracers
1634                     do k=1,levs
1635                       do i=1,im
1636                         clw(i,k,2+n) = gq0(i,k,n+trc_shft)
1637                       enddo
1638                     enddo
1639                 enddo
1640               endif
1641     
1642             endif
1643           endif   ! end if_ras
1644     
1645           do i = 1, im
1646             ktop(i)  = 1
1647             kbot(i)  = levs
1648           enddo
1649     
1650     !  --- ...  calling precipitation processes
1651     
1652     !     do i = 1, im
1653     !       work1(i) = (log(coslat(i) / (nlons(i)*latg)) - dxmin) * dxinv
1654     !       work1(i) = max(0.0, min(1.0,work1(i)))
1655     !       work2(i) = 1.0 - work1(i)
1656     !     enddo
1657     
1658     !  --- ...  calling convective parameterization
1659     
1660           if (ntcw > 0) then
1661     
1662             do k = 1, levs
1663               do i = 1, im
1664                 rhc(i,k) = rhbbot - (rhbbot-rhbtop) * (1.0-prslk(i,k))
1665                 rhc(i,k) = rhc_max * work1(i) + rhc(i,k) * work2(i)
1666                 rhc(i,k) = max(0.0, min(1.0,rhc(i,k)))
1667               enddo
1668             enddo
1669     
1670             if (num_p3d == 3) then    ! call brad ferrier's microphysics
1671               do i = 1, im
1672                 flgmin_l(i) = flgmin(1)*work1(i) + flgmin(2)*work2(i)
1673               enddo
1674     
1675     !  --- ...  algorithm to separate different hydrometeor species
1676     
1677               do k = 1, levs
1678                 do i = 1, im
1679                   wc     = gq0(i,k,ntcw)
1680                   qi     = 0.
1681                   qr     = 0.
1682                   qw     = 0.
1683                   f_ice  = max(0.0, min(1.0, phy_f3d(i,k,1)))
1684                   f_rain = max(0.0, min(1.0, phy_f3d(i,k,2)))
1685     
1686                   qi = f_ice*wc
1687                   qw = wc-qi
1688                   if (qw > 0.0) then
1689                     qr = f_rain*qw
1690                     qw = qw-qr
1691                   endif
1692     
1693     !             if (f_ice >= 1.0) then
1694     !               qi = wc
1695     !             elseif (f_ice <= 0.0) then
1696     !               qw = wc
1697     !             else
1698     !               qi = f_ice*wc
1699     !               qw = wc-qi
1700     !             endif
1701     
1702     !             if (qw > 0.0 .and. f_rain > 0.0) then
1703     !               if (f_rain >= 1.0) then
1704     !                 qr = qw
1705     !                 qw = 0.0
1706     !               else
1707     !                 qr = f_rain*qw
1708     !                 qw = qw-qr
1709     !               endif
1710     !             endif
1711     
1712                   qr_col(i,k) = qr
1713     !             clw(i,k)    = qi + qw
1714                   clw(i,k,1)  = qi
1715                   clw(i,k,2)  = qw
1716     
1717     !  --- ...  array to track fraction of "cloud" in the form of ice
1718     
1719     !             if (qi+qw > epsq) then
1720     !               fc_ice(i,k) = qi / (qi+qw)
1721     !             else
1722     !               fc_ice(i,k) = 0.0
1723     !             endif
1724                 enddo
1725               enddo
1726     
1727             else   ! if_num_p3d
1728     
1729               do i = 1, im
1730                 psautco_l(i) = psautco(1)*work1(i) + psautco(2)*work2(i)
1731                 prautco_l(i) = prautco(1)*work1(i) + prautco(2)*work2(i)
1732               enddo
1733               do k = 1, levs
1734                 do i = 1, im
1735                   clw(i,k,1) = gq0(i,k,ntcw)
1736                 enddo
1737               enddo
1738     
1739             endif  ! end if_num_p3d
1740     
1741           else    ! if_ntcw
1742     
1743             do i = 1, im
1744               psautco_l(i) = psautco(1)*work1(i) + psautco(2)*work2(i)
1745               prautco_l(i) = prautco(1)*work1(i) + prautco(2)*work2(i)
1746             enddo
1747             do k = 1, levs
1748               do i = 1, im
1749                 rhc(i,k) = 1.0
1750               enddo
1751             enddo
1752     
1753           endif   ! end if_ntcw
1754     !
1755           if (.not. ras) then
1756     
1757             if (.not. newsas) then
1758     
1759               call sascnv(im,ix,levs,jcap,dtp,del,prsl,pgr,phil,            &
1760          &                clw,gq0,gt0,gu0,gv0,cld1d,                        &
1761          &                rain1,kbot,ktop,kcnv,slmsk,                       &
1762          &                vvel,rann,ncld,ud_mf,dd_mf,dt_mf)
1763     !hchuang code change 03/03/08 [r1L] add SAS modification of mass flux
1764     !    &                vvel,rann,ncld)
1765     
1766             else
1767     
1768               call sascnvn(im,ix,levs,jcap,dtp,del,prsl,pgr,phil,           &
1769          &                clw,gq0,gt0,gu0,gv0,cld1d,                        &
1770          &                rain1,kbot,ktop,kcnv,slmsk,                       &
1771          &                VVEL,ncld,ud_mf,dd_mf,dt_mf)
1772     !hchuang code change 03/03/08 [r1L] add SAS modification of mass flux
1773     !    &                vvel,rann,ncld)
1774     
1775             endif
1776     
1777     !       if (lprnt) print *,' rain1=',rain1(ipr),' rann=',rann(ipr,1)
1778     
1779           else    ! if_not_newsas
1780     
1781     !       if  (lprnt) print *,' calling ras for kdt=',kdt,' me=',me       &
1782     !    &,                     ' lprnt=',lprnt
1783     
1784             if (ccwf(1) >= 0.0 .or. ccwf(2) >= 0 ) then
1785               do i=1,im
1786                 ccwfac(i) = ccwf(1)*work1(i) + ccwf(2)*work2(i)
1787                 dlqfac(i) = dlqf(1)*work1(i) + dlqf(2)*work2(i)
1788               enddo
1789             else
1790               ccwfac = -999.0
1791               dlqfac = 0.0
1792             endif
1793     
1794             call rascnv(im,     ix,    levs,   dtp, dtf, rann               &
1795          &,             gt0,    gq0,   gu0,    gv0, clw, tottracer          &
1796          &,             prsi,   prsl,   prsik,  prslk, phil,  phii          &
1797          &,             kpbl,   cd,     rain1,  kbot,  ktop,  kcnv          &
1798          &,             phy_f2d(1,num_p2d), flipv, pa2mb                    &
1799          &,             me, garea, lmh, ccwfac, nrcm, rhc                   &
1800          &,             ud_mf, dd_mf, dt_mf, dlqfac, lprnt, ipr, kdt, fscav)
1801     
1802     !  --- ...  check print
1803     
1804     !       if (lprnt) print *,' rain1=',rain1(ipr),' rann=',rann(ipr,1)
1805     !    &,' lat=',lat
1806     !       do i = 1, im
1807     !         if (tsea(i) > 380.0 .or. tsea(i) < 10) then
1808     !           print *,' tsea=', tsea(i),' i=',i,' lat=',lat,              &
1809     !    &        ' kdt=',kdt,' xlon=',xlon(i),' xlat=',xlat(i),' slmsk=',  &
1810     !    &        slmsk(i),' me=',me
1811     !           stop
1812     !         endif
1813     !       enddo
1814     
1815     !       do k = 1, levs
1816     !         do i = 1, im
1817     !           if (gt0(i,k) > 330.0 .or. gt0(i,k) < 80.0) then             &
1818     !             print *,' gt0=', gt0(i,k),' i=',i,' k=',k,' lat=',lat,    &
1819     !    &          ' kdt=',kdt,' xlon=',xlon(i),' xlat=',xlat(i)
1820     !             stop
1821     !           endif
1822     
1823     !           if (gq0(i,k,1) > 1.0 ) then
1824     !             print *,' gq0=', gq0(i,k,1),' i=',i,' k=',k,' lat=',lat,  &
1825     !    &          ' kdt=',kdt
1826     !             stop
1827     !           endif
1828     !         enddo
1829     !       enddo
1830     
1831     !       if (lprnt) print *,' returning from ras for kdt=', kdt,         &
1832     !    &                     ' me=',me,' lat=',lat
1833     
1834     !  --- ...  end check print
1835     
1836             cld1d = 0
1837     
1838     !  --- ...  update the tracers due to convective transport
1839     
1840             if (tottracer > 0) then
1841               if (ntoz > 0) then                         ! for ozone
1842                 do k=1,levs
1843                   do i=1,im
1844                     gq0(i,k,ntoz) = clw(i,k,3)
1845                   enddo
1846                 enddo
1847     
1848                 if (tracers > 0) then                    ! for other tracers
1849                   do n = 1, tracers
1850                     do k=1,levs
1851                       do i=1,im
1852                         gq0(i,k,n+trc_shft) = clw(i,k,3+n)
1853                       enddo
1854                     enddo
1855                   enddo
1856                 endif
1857               else
1858                 do n = 1, tracers
1859                   do k=1,levs
1860                     do i=1,im
1861                       gq0(i,k,n+trc_shft) = clw(i,k,2+n)
1862                     enddo
1863                   enddo
1864                 enddo
1865               endif
1866             endif
1867           endif   ! end if_not_ras
1868     !
1869           if (lssav) then
1870             do i = 1, im
1871               cldwrk(i) = cldwrk(i) + cld1d(i) * dtf
1872             enddo
1873     
1874             if (ldiag3d) then
1875               do k = 1, levs
1876                 do i = 1, im
1877                   dt3dt(i,k,4) = dt3dt(i,k,4) + (gt0(i,k)-dtdt(i,k)) * frain
1878     !             dq3dt(i,k,2) = dq3dt(i,k,2) + (gq0(i,k,1)-dqdt(i,k,1))    &
1879     !    &                                                           * frain
1880                   du3dt(i,k,3) = du3dt(i,k,3) + (gu0(i,k)-dudt(i,k)) * frain
1881                   dv3dt(i,k,3) = dv3dt(i,k,3) + (gv0(i,k)-dvdt(i,k)) * frain
1882     
1883                 enddo
1884               enddo
1885             endif
1886     ! update dqdt_v to include moisture tendency due to deep convection
1887             if (lgocart) then
1888               do k = 1, levs
1889                 do i = 1, im
1890                   tem          = (gq0(i,k,1)-dqdt(i,k,1)) * frain
1891                   dqdt_v(i,k)  = dqdt_v(i,k)  + tem
1892                 enddo
1893               enddo
1894             endif
1895             if (ldiag3d .or. lggfs3d) then
1896               do k = 1, levs
1897                 do i = 1, im
1898                   tem          = (gq0(i,k,1)-dqdt(i,k,1)) * frain
1899                   dq3dt(i,k,2) = dq3dt(i,k,2) + tem
1900     !             dqdt_v(i,k)  = dqdt_v(i,k)  + tem
1901                   upd_mf(i,k)  = upd_mf(i,k)  + ud_mf(i,k) * (con_g*frain)
1902                   dwn_mf(i,k)  = dwn_mf(i,k)  + dd_mf(i,k) * (con_g*frain)
1903                   det_mf(i,k)  = det_mf(i,k)  + dt_mf(i,k) * (con_g*frain)
1904                 enddo
1905               enddo
1906             endif
1907           endif   ! end if_lssav
1908     
1909           do i = 1, im
1910             rainc(i) = frain * rain1(i)
1911           enddo
1912     
1913           if (lssav) then
1914             do i = 1, im
1915               cnvprcp(i) = cnvprcp(i) + rainc(i)
1916             enddo
1917           endif
1918     !
1919           if (cnvgwd) then         !        call convective gravity wave drag
1920     
1921     !  --- ...  calculate maximum convective heating rate            qmax [k/s]
1922     !           cuhr = temperature change due to deep convection
1923     
1924             do i = 1, im
1925               qmax(i)   = 0.
1926               cumabs(i) = 0.
1927             enddo
1928     
1929             do k = 1, levs
1930               do i = 1, im
1931     !           cuhr(i,k) = (gt0(i,k)-dtdt(i,k)) / dtf
1932                 cuhr(i,k) = (gt0(i,k)-dtdt(i,k)) / dtp    ! moorthi
1933     
1934                 cumchr(i,k) = 0.
1935                 gwdcu(i,k)  = 0.
1936                 gwdcv(i,k)  = 0.
1937                 diagn1(i,k) = 0.
1938                 diagn2(i,k) = 0.
1939     
1940                 if (k >= kbot(i) .and. k <= ktop(i)) then
1941                   qmax(i)     = max(qmax(i),cuhr(i,k))
1942                   cumabs(i)   = cuhr(i,k) + cumabs(i)
1943                 endif
1944               enddo
1945             enddo
1946     
1947             do i = 1, im
1948               do k = kbot(i), ktop(i)
1949                 do k1 = kbot(i), k
1950                   cumchr(i,k) = cuhr(i,k1) + cumchr(i,k)
1951                 enddo
1952     
1953                 cumchr(i,k) = cumchr(i,k) / cumabs(i)
1954               enddo
1955             enddo
1956     
1957     !  --- ...  check print
1958     
1959             if (lprnt) then
1960               if (kbot(ipr) <= ktop(ipr)) then
1961                 write(*,*) 'kbot <= ktop     for (lat,lon) = ',             &
1962          &            xlon(ipr)*57.29578,xlat(ipr)*57.29578
1963                 write(*,*) 'kcnv kbot ktop qmax dlength  ',kcnv(ipr),       &
1964          &            kbot(ipr),ktop(ipr),(86400.*qmax(ipr)),dlength(ipr)
1965                 write(*,9000) kdt
1966      9000       format(/,3x,'k',5x,'cuhr(k)',4x,'cumchr(k)',5x,             &
1967          &            'at kdt = ',i4,/)
1968     
1969                 do k = ktop(ipr), kbot(ipr),-1
1970                   write(*,9010) k,(86400.*cuhr(ipr,k)),(100.*cumchr(ipr,k))
1971      9010         format(2x,i2,2x,f8.2,5x,f6.0)
1972                 enddo
1973               endif
1974     
1975     !         print *,' before gwdc in gbphys fhour ',fhour
1976     
1977               if (fhour >= fhourpr) then
1978                 print *,' before gwdc in gbphys start print'
1979                 write(*,*) 'fhour ix im levs = ',fhour,ix,im,levs
1980                 print *,'dtp  dtf  = ',dtp,dtf
1981     
1982                 write(*,9100)
1983      9100       format(//,14x,'pressure levels',//                          &
1984          &             ' ilev',7x,'prsi',8x,'prsl',8x,'delp',/)
1985     
1986                 k = levs + 1
1987                 write(*,9110) k,(10.*prsi(ipr,k))
1988      9110       format(i4,2x,f10.3)
1989     
1990                 do k = levs, 1, -1
1991                   write(*,9120) k,(10.*prsl(ipr,k)),(10.*del(ipr,k))
1992                   write(*,9110) k,(10.*prsi(ipr,k))
1993                 enddo
1994      9120       format(i4,12x,2(2x,f10.3))
1995     
1996                 write(*,9130)
1997      9130       format(//,10x,'before gwdc in gbphys',//,' ilev',6x,        &
1998          &             'ugrs',9x,'gu0',8x,'vgrs',9x,'gv0',8x,               &
1999          &             'tgrs',9x,'gt0',8x,'gt0b',8x,'dudt',8x,'dvdt',/)
2000     
2001                 do k = levs, 1, -1
2002                   write(*,9140) k,ugrs(ipr,k),gu0(ipr,k),                   &
2003          &                        vgrs(ipr,k),gv0(ipr,k),                   &
2004          &                        tgrs(ipr,k),gt0(ipr,k),dtdt(ipr,k),       &
2005          &                        dudt(ipr,k),dvdt(ipr,k)
2006                 enddo
2007      9140       format(i4,9(2x,f10.3))
2008     
2009                 print *,' before gwdc in gbphys end print'
2010               endif
2011             endif   ! end if_lprnt
2012     
2013     !  --- ...  end check print
2014     
2015             call gwdc(im, ix, im, levs, lat, ugrs, vgrs, tgrs, qgrs,        &
2016          &           prsl, prsi, del, qmax, cumchr, ktop, kbot, kcnv,       &
2017          &          gwdcu, gwdcv,con_g,con_cp,con_rd,con_fvirt, dlength,    &
2018          &          lprnt, ipr, fhour,                                      &
2019          &          dusfcg,dvsfcg,diagn1,diagn2)
2020     
2021             if (lprnt) then
2022               if (fhour >= fhourpr) then
2023                 print *,' after gwdc in gbphys start print'
2024     
2025                 write(*,9131)
2026      9131       format(//,10x,'after gwdc in gbphys',//,' ilev',6x,         &
2027          &             'ugrs',9x,'gu0',8x,'vgrs',9x,'gv0',8x,               &
2028          &             'tgrs',9x,'gt0',8x,'gt0b',7x,'gwdcu',7x,'gwdcv',/)
2029     
2030                 do k = levs, 1, -1
2031                   write(*,9141) k,ugrs(ipr,k),gu0(ipr,k),                   &
2032          &                        vgrs(ipr,k),gv0(ipr,k),                   &
2033          &                        tgrs(ipr,k),gt0(ipr,k),dtdt(ipr,k),       &
2034          &                        gwdcu(ipr,k),gwdcv(ipr,k)
2035                 enddo
2036      9141       format(i4,9(2x,f10.3))
2037     
2038                 print *,' after gwdc in gbphys end print'
2039               endif
2040             endif
2041     
2042     !  --- ...  write out cloud top stress and wind tendencies
2043     
2044             if (lssav) then
2045               do i = 1, im
2046                 dugwd(i) = dugwd(i) + dusfcg(i)*dtf
2047                 dvgwd(i) = dvgwd(i) + dvsfcg(i)*dtf
2048               enddo
2049     
2050               if (ldiag3d) then
2051                 do k = 1, levs
2052                   do i = 1, im
2053                     du3dt(i,k,4) = du3dt(i,k,4) + gwdcu(i,k)  * dtf
2054                     dv3dt(i,k,4) = dv3dt(i,k,4) + gwdcv(i,k)  * dtf
2055     !               du3dt(i,k,2) = du3dt(i,k,2) + diagn1(i,k) * dtf
2056     !               dv3dt(i,k,2) = dv3dt(i,k,2) + diagn2(i,k) * dtf
2057                   enddo
2058                 enddo
2059               endif
2060             endif   ! end if_lssav
2061     
2062     !  --- ...  update the wind components with  gwdc tendencies
2063     
2064             do k = 1, levs
2065               do i = 1, im
2066                 gu0(i,k)   = gu0(i,k)   + gwdcu(i,k)   * dtp
2067                 gv0(i,k)   = gv0(i,k)   + gwdcv(i,k)   * dtp
2068               enddo
2069             enddo
2070     
2071             if (lprnt) then
2072               if (fhour >= fhourpr) then
2073                 print *,' after tendency gwdc in gbphys start print'
2074     
2075                 write(*,9132)
2076      9132       format(//,10x,'after tendency gwdc in gbphys',//,' ilev',6x,&
2077          &             'ugrs',9x,'gu0',8x,'vgrs',9x,'gv0',8x,               &
2078          &             'tgrs',9x,'gt0',8x,'gt0b',7x,'gwdcu',7x,'gwdcv',/)
2079     
2080                 do k = levs, 1, -1
2081                   write(*,9142) k,ugrs(ipr,k),gu0(ipr,k),vgrs(ipr,k),       &
2082          &              gv0(ipr,k),tgrs(ipr,k),gt0(ipr,k),dtdt(ipr,k),      &
2083          &              gwdcu(ipr,k),gwdcv(ipr,k)
2084                 enddo
2085      9142       format(i4,9(2x,f10.3))
2086     
2087                 print *,' after tendency gwdc in gbphys end print'
2088               endif
2089             endif
2090     
2091           endif   ! end if_cnvgwd (convective gravity wave drag)
2092     
2093           if (ldiag3d) then
2094             do k = 1, levs
2095               do i = 1, im
2096                 dtdt(i,k)   = gt0(i,k)
2097     !           dqdt(i,k,1) = gq0(i,k,1)
2098               enddo
2099             enddo
2100           endif
2101           if (ldiag3d .or. lggfs3d) then
2102             do k = 1, levs
2103               do i = 1, im
2104                 dqdt(i,k,1) = gq0(i,k,1)
2105               enddo
2106             enddo
2107           endif
2108     
2109           if (shal_cnv) then
2110             if (.not. sashal) then
2111     
2112     !         if (lprnt) print *,' levshcm=',levshcm,' gt0sh=',gt0(ipr,:)
2113     !    &,  ' lat=',lat
2114     
2115               if (.not. mstrat) then
2116                 call shalcvt3(im,ix,levshcm,dtp,del,prsi,prsl,prslk,        &
2117          &                    kcnv,gq0,gt0)                             ! for pry
2118               else
2119                 call shalcv(im,ix,levshcm,dtp,del,prsi,prsl,prslk,kcnv,     &
2120          &                  gq0,gt0,levshc,phil,kinver,ctei_r,ctei_rml      &
2121          &,                                                lprnt,ipr)
2122     
2123     !           if (lprnt) print *,' levshcm=',levshcm,' gt0sha=',gt0(ipr,:)
2124               endif
2125     
2126             else
2127     
2128               call shalcnv(im,ix,levs,jcap,dtp,del,prsl,pgr,phil,           &
2129          &                 clw,gq0,gt0,gu0,gv0,                             &
2130          &                 rain1,kbot,ktop,kcnv,slmsk,                      &
2131          &                 vvel,ncld,hpbl,hflx,evap,ud_mf,dt_mf)
2132     
2133               do i = 1, im
2134                 raincs(i) = frain * rain1(i)
2135               enddo
2136     
2137               if (lssav) then
2138                 do i = 1, im
2139                   cnvprcp(i) = cnvprcp(i) + raincs(i)
2140                 enddo
2141               endif
2142     
2143               do i = 1, im
2144                 rainc(i) = rainc(i) + raincs(i)
2145               enddo
2146     
2147               if (lssav) then
2148                 do i = 1, im
2149                   cldwrk(i) = cldwrk(i) + cld1d(i) * dtf
2150                 enddo
2151               endif
2152     
2153             endif   ! end if_not_sashal
2154           endif   ! end if_shal_cnv
2155     
2156           if (lssav) then
2157     ! update dqdt_v to include moisture tendency due to shallow convection
2158             if (lgocart) then
2159               do k = 1, levs
2160                 do i = 1, im
2161                   tem          = (gq0(i,k,1)-dqdt(i,k,1)) * frain
2162                   dqdt_v(i,k)  = dqdt_v(i,k)  + tem
2163                 enddo
2164               enddo
2165             endif
2166             if (ldiag3d) then
2167               do k = 1, levs
2168                 do i = 1, im
2169                   dt3dt(i,k,5) = dt3dt(i,k,5) + (gt0(i,k)-dtdt(i,k)) * frain
2170     !             dq3dt(i,k,3) = dq3dt(i,k,3) + (gq0(i,k,1)-dqdt(i,k,1))    &
2171     !    &                                                           * frain
2172                 enddo
2173               enddo
2174               do k = 1, levs
2175                 do i = 1, im
2176                   dtdt(i,k)   = gt0(i,k)
2177     !             dqdt(i,k,1) = gq0(i,k,1)
2178                 enddo
2179               enddo
2180             endif
2181             if (ldiag3d .or. lggfs3d) then
2182               do k = 1, levs
2183                 do i = 1, im
2184                   tem          = (gq0(i,k,1)-dqdt(i,k,1)) * frain
2185                   dq3dt(i,k,3) = dq3dt(i,k,3) + tem
2186     !             dqdt_v(i,k)  = dqdt_v(i,k)  + tem
2187     !             dqdt_v(i,k)  = dqdt_v(i,k)  * (1000.0/dtf)
2188                 enddo
2189               enddo
2190               do k = 1, levs
2191                 do i = 1, im
2192                   dqdt(i,k,1) = gq0(i,k,1)
2193                 enddo
2194               enddo
2195             endif
2196           endif   ! end if_lssav
2197     
2198           do k = 1, levs
2199             do i = 1, im
2200               if (clw(i,k,2) <= -999.0) clw(i,k,2) = 0.0
2201             enddo
2202           enddo
2203     
2204           if (ntcw > 0) then
2205     
2206             if (num_p3d == 3) then    ! call brad ferrier's microphysics
2207     
2208               do k = 1, levs
2209                 do i = 1, im
2210     
2211     !             qi = clw(i,k)*fc_ice(i,k)
2212     !             qw = clw(i,k) - qi
2213                   qi = clw(i,k,1)
2214                   qw = clw(i,k,2)
2215     
2216     !  --- ...  algorithm to combine different hydrometeor species
2217     
2218     !             gq0(i,k,ntcw) = max(epsq, qi+qw+qr_col(i,k))
2219                   gq0(i,k,ntcw) = qi + qw + qr_col(i,k)
2220     
2221                   if (qi <= epsq) then
2222                     phy_f3d(i,k,1) = 0.
2223                   else
2224                     phy_f3d(i,k,1) = qi/gq0(i,k,ntcw)
2225                   endif
2226     
2227                   if (qr_col(i,k) <= epsq) then
2228                     phy_f3d(i,k,2) = 0.
2229                   else
2230                     phy_f3d(i,k,2) = qr_col(i,k) / (qw+qr_col(i,k))
2231                   endif
2232     
2233                 enddo
2234               enddo
2235     
2236             else    ! if_num_p3d
2237     
2238               do k = 1, levs
2239                 do i = 1, im
2240     !             gq0(i,k,ntcw) = clw(i,k)   + gq0(i,k,ntcw)
2241                   gq0(i,k,ntcw) = clw(i,k,1) + clw(i,k,2)
2242     !             gq0(i,k,ntcw) = clw(i,k,1)               ! for pry
2243                 enddo
2244               enddo
2245     
2246             endif   ! end if_num_p3d
2247     
2248           else    ! if_ntcw
2249     
2250             do k = 1, levs
2251               do i = 1, im
2252                 clw(i,k,1) = clw(i,k,1) + clw(i,k,2)
2253               enddo
2254             enddo
2255     
2256           endif   ! end if_ntcw
2257     
2258           call cnvc90(clstp, im,   ix,   rainc, kbot, ktop, levs, prsi,     &
2259          &            acv,   acvb, acvt, cv,    cvb,  cvt)
2260     
2261           if (ncld == 0) then
2262     
2263             call lrgscl(ix,im,levs,dtp,gt0,gq0,prsl,del,prslk,rain1,clw)
2264     
2265           elseif (ncld == 1) then
2266             if (moist_adj) then
2267     !
2268     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2269     !         if (me .eq. 0) then
2270     !           sumq = 0.0
2271     !           DO K=1,LEVS
2272     !             do i=1,im
2273     !               sumq = sumq - (gq0(i,k,1)+gq0(i,k,ntcw)) * del(i,k)
2274     !             enddo
2275     !           enddo
2276     !         endif
2277     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2278     !
2279     !         To call moist convective adjustment
2280     !
2281     !         if (lprnt) then
2282     !           print *,' prsl=',prsl(ipr,:)
2283     !           print *,' del=',del(ipr,:)
2284     !           print *,' gt0b=',gt0(ipr,:)
2285     !           print *,' gq0b=',gq0(ipr,:,1)
2286     !         endif
2287     
2288               call mstcnv(im,ix,levs,dtp,gt0,gq0,prsl,del,prslk,rain1
2289          &,                            gq0(1,1,ntcw), rhc, lprnt,ipr)
2290     
2291     !         if (lprnt) then
2292     !           print *,' rain1=',rain1(ipr),' rainc=',rainc(ipr)
2293     !           print *,' gt0a=',gt0(ipr,:)
2294     !           print *,' gq0a=',gq0(ipr,:,1)
2295     !         endif
2296               do i=1,im
2297                 rainc(i) = rainc(i) + frain * rain1(i)
2298               enddo
2299               if(lssav) then
2300                 do i=1,im
2301                   cnvprcp(i) = cnvprcp(i) + rain1(i) * frain
2302                 enddo
2303     ! update dqdt_v to include moisture tendency due to surface processes
2304     ! dqdt_v : instaneous moisture tendency (kg/kg/sec)
2305                 if (lgocart) then
2306                   do k=1,levs
2307                     do i=1,im
2308                       tem = (gq0(i,k,1)-dqdt(i,k,1)) * frain
2309                       dqdt_v(i,k) = dqdt_v(i,k) + tem
2310                       dqdt_v(i,k) = dqdt_v(i,k) / dtf
2311                     enddo
2312                   enddo
2313                 endif
2314                 if (ldiag3d) then
2315                   do k=1,levs
2316                     do i=1,im
2317                       dt3dt(i,k,4) = dt3dt(i,k,4) + (gt0(i,k)-dtdt(i,k))
2318          &                                        * frain
2319     !                 dq3dt(i,k,2) = dq3dt(i,k,2) + (gq0(i,k,1)-dqdt(i,k,1))
2320     !    &                                        * frain
2321                     enddo
2322                   enddo
2323                 endif
2324                 if (ldiag3d .or. lggfs3d) then
2325                   do k=1,levs
2326                     do i=1,im
2327                       dq3dt(i,k,2) = dq3dt(i,k,2) + (gq0(i,k,1)-dqdt(i,k,1))
2328          &                                        * frain
2329                     enddo
2330                   enddo
2331                 endif
2332                endif
2333     !
2334               if (ldiag3d) then
2335                 do k=1,levs
2336                   do i=1,im
2337                     dtdt(i,k)   = gt0(i,k)
2338     !               dqdt(i,k,1) = gq0(i,k,1)
2339                   enddo
2340                 enddo
2341               endif
2342               if (ldiag3d .or. lggfs3d) then
2343                 do k=1,levs
2344                   do i=1,im
2345                     dqdt(i,k,1) = gq0(i,k,1)
2346                   enddo
2347                 enddo
2348               ENDIF
2349     !
2350     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2351     !         if (me .eq. 0) then
2352     !           DO K=1,LEVS
2353     !             do i=1,im
2354     !               sumq = sumq + (gq0(i,k,1)+gq0(i,k,ntcw)) * del(i,k)
2355     !             enddo
2356     !           enddo
2357     !           sumr = 0.0
2358     !           do i=1,im
2359     !             sumr = sumr + rain1(i)
2360     !           enddo
2361     !           sumq = sumq * 1000.0 / grav
2362     !           sumr = sumr *1000
2363     !           print *,' after MCN: sumq=',sumq,' sumr=',sumr, ' kdt=',kdt
2364     !         endif
2365     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2366     !
2367             endif               !       moist convective adjustment over
2368     !
2369     
2370             if (num_p3d == 3) then    ! call brad ferrier's microphysics
2371     
2372               do i = 1, im
2373                 xncw(i) = ncw(2) * work1(i) + ncw(1) * work2(i)
2374               enddo
2375     
2376               if (kdt == 1 .and. abs(xlon(1)) < 0.0001) then
2377                 print *,' xncw=',xncw(1),' rhc=',rhc(1,1),' work1=',work1(1)&
2378          &,         ' work2=',work2(1),' flgmin=',flgmin_l(1)               &
2379          &,         ' lon=',xlon(1) * 57.29578,' lat=',lat,' me=',me
2380     !    &,         ' lon=',xlon(1) * 57.29578,' lat=',xlat(1) * 57.29578
2381     !    &,         ' kinver=',kinver(1)
2382               endif
2383     
2384     !         if (lprnt) print *,' ipr=',ipr,' gt0_gsmb=',gt0(ipr,:)        &
2385     !    &,         ' xlon=',xlon(ipr),' xlat=',xlat(ipr)
2386     
2387               call gsmdrive(im, ix, levs, dtp, prsl, del                    &
2388          &,                 gt0, gq0(1,1,1), gq0(1,1,ntcw), slmsk           &
2389          &,                 phy_f3d(1,1,1),  phy_f3d(1,1,2)                 &
2390          &,                 phy_f3d(1,1,3), rain1, sr, con_g                &
2391          &,                 con_hvap, hsub,con_cp, rhc, xncw, flgmin_l      &
2392          &,                 me, lprnt, ipr)
2393     
2394     !         if (lprnt) print *,' ipr=',ipr,' gt0_gsma=',gt0(ipr,:)
2395     
2396             elseif (num_p3d == 4) then ! call zhao/carr/sundqvist microphysics
2397     
2398               call gscond(im, ix, levs, dtp, prsl, pgr,                     &
2399          &                gq0(1,1,1), gq0(1,1,ntcw), gt0,                   &
2400          &                phy_f3d(1,1,1), phy_f3d(1,1,2), phy_f2d(1,1),     &
2401          &                phy_f3d(1,1,3), phy_f3d(1,1,4), phy_f2d(1,2),     &
2402          &                rhc,lprnt, ipr)
2403     
2404               call precpd(im, ix, levs, dtp, del, prsl, pgr,                &
2405          &                gq0(1,1,1), gq0(1,1,ntcw), gt0, rain1,            &
2406          &                rainp, rhc, psautco_l, prautco_l, evpco,          &
2407          &                lprnt, ipr)
2408     
2409             endif   ! end if_num_p3d
2410     
2411           endif   ! end if_ncld
2412     
2413     !     if (lprnt) print *,' rain1=',rain1(ipr),' rainc=',rainc(ipr)
2414     
2415     
2416           do i = 1, im
2417             rainl(i) = frain    * rain1(i)
2418             rain(i)  = rainc(i) + rainl(i)
2419           enddo
2420     
2421     !  --- ...  coupling insertion
2422     
2423           if (lssav_cc) then
2424             precr_cc(1:im) = precr_cc(1:im) + rain(1:im)
2425           endif
2426     
2427     !  --- ...  end coupling insertion
2428     
2429           if (cal_pre) then
2430     !       HCHUANG: add dominant precipitation type algorithm
2431     
2432             call calpreciptype(kdt,nrcm,im,ix,LEVS,LEVS+1,rann,
2433          &    xlat,xlon,gt0,gq0,prsl,prsi,RAIN,
2434          &    phii,num_p3d,tsea,sr,phy_f3d(1,1,3),    ! input
2435          &    DOMR,DOMZR,DOMIP,DOMS)                  ! output
2436     
2437     !
2438     !        if (lprnt) print*,'debug calpreciptype: DOMR,DOMZR,DOMIP,DOMS '
2439     !     &,DOMR(ipr),DOMZR(ipr),DOMIP(ipr),DOMS(ipr)
2440     !        do i=1,im
2441     !         if (abs(xlon(i)*57.29578-114.0) .lt. 0.2  .and.
2442     !     &    abs(xlat(i)*57.29578-40.0) .lt. 0.2)
2443     !     &    print*,'debug calpreciptype: DOMR,DOMZR,DOMIP,DOMS ',
2444     !     &    DOMR(i),DOMZR(i),DOMIP(i),DOMS(i)
2445     !       end do
2446     !       HCHUANG: use new precipitation type to decide snow flag for LSM snow accumulation
2447             do i=1,im
2448               if(DOMS(i) >0.0 .or. DOMIP(i)>0.0)then
2449                 SRFLAG(i) = 1.
2450               else
2451                 SRFLAG(i) = 0.
2452               end if
2453             enddo
2454           endif
2455     
2456           if (lssav) then
2457             do i = 1, im
2458               totprcp(i) = totprcp(i) + rain(i)
2459             enddo
2460     
2461             if (ldiag3d) then
2462               do k = 1, levs
2463                 do i = 1, im
2464                   dt3dt(i,k,6) = dt3dt(i,k,6) + (gt0(i,k)-dtdt(i,k)) * frain
2465     !             dq3dt(i,k,4) = dq3dt(i,k,4) + (gq0(i,k,1)-dqdt(i,k,1))    &
2466     !    &                                                           * frain
2467                 enddo
2468               enddo
2469             endif
2470             if (ldiag3d .or. lggfs3d) then
2471               do k = 1, levs
2472                 do i = 1, im
2473                   dq3dt(i,k,4) = dq3dt(i,k,4) + (gq0(i,k,1)-dqdt(i,k,1))    &
2474          &                                                           * frain
2475                 enddo
2476               enddo
2477             endif
2478     !hchuang code change [+8L] 10/08/2008 rnp is 1.E+3 m/dtp, with mutiply by frain dtf/dtp
2479     !              rnp becomes rain (km) amount when output it will become
2480     !              rnp/(total time)  that is km/s
2481             if (lggfs3d) then
2482               do k = 1, levs
2483                 do i = 1, im
2484                   rnp(i,k) = rainp(i,k) * frain ! 3D large scale liqid rain amount
2485                                                 ! of one DTF time period in mm
2486                 enddo
2487               enddo
2488             endif
2489           endif
2490     
2491     !  --- ...  estimate t850 for rain-snow decision
2492     
2493           do i = 1, im
2494             t850(i) = gt0(i,1)
2495           enddo
2496     
2497           do k = 1, levs - 1
2498             do i = 1, im
2499               if (prsl(i,k) > p850 .and. prsl(i,k+1) <= p850) then
2500                 t850(i) = gt0(i,k) - (prsl(i,k)-p850)                       &
2501          &              / (prsl(i,k)-prsl(i,k+1)) * (gt0(i,k)-gt0(i,k+1))
2502               endif
2503             enddo
2504           enddo
2505     
2506     !  --- ...  lu: snow-rain detection is performed in land/sice module
2507     
2508           if (cal_pre) then ! HCHUANG: new precip type algorithm defines srflag
2509             do i = 1, im
2510               tprcp(i) = rain(i)            ! clu: rain -> tprcp
2511               if (t850(i) <= 273.16) then
2512     !  --- ...  wei: when call osu lsm, neutral impact, for code consistency
2513                 if (lsm == 0 .and. slmsk(i) /= 0.0) then
2514                   weasd(i)  = weasd(i) + 1.e3*rain(i)
2515                   tprcp(i)  = 0.
2516                 endif
2517               endif
2518             enddo
2519           else
2520             do i = 1, im
2521               tprcp(i) = rain(i)            ! clu: rain -> tprcp
2522               srflag(i) = 0.                ! clu: default srflag as 'rain' (i.e. 0)
2523     
2524               if (t850(i) <= 273.16) then
2525                 srflag(i) = 1.              ! clu: set srflag to 'snow' (i.e. 1)
2526     !  --- ...  wei: when call osu lsm, neutral impact, for code consistency
2527                 if (lsm == 0 .and. slmsk(i) /= 0.0) then
2528                   weasd(i)  = weasd(i) + 1.e3*rain(i)
2529                   tprcp(i)  = 0.
2530                 endif
2531               endif
2532             enddo
2533           endif
2534     
2535     !     if (lprnt) print *,' tprcp=',tprcp(ipr),' rain=',rain(ipr)
2536     
2537     !  --- ...  cpl insertion
2538     
2539           do i = 1, im
2540     !       if (t850(i) <= 273.16 .and. slmsk(i) /= 0.) then
2541             if (t850(i) <= 273.16) then
2542                lprec_cc(i) = 0.0
2543                snw_cc(i)   = rain(i)
2544             else
2545                lprec_cc(i) = rain(i)
2546                snw_cc(i)   = 0.0
2547             endif
2548           enddo
2549     
2550     !  --- ...  wei: when call osu lsm
2551     !           update soil moisture and canopy water after precipitation computaion
2552     
2553           if (lsm == 0) then
2554     
2555             call progt2(im,lsoil,rhscnpy,rhsmc,ai,bi,cci,smsoil,            &
2556          &              slmsk,canopy,tprcp,runof,snowmt,                    &
2557          &              zsoil,soiltyp,sigmaf,dtf,me)
2558     
2559     !  --- ...  wei: let soil liquid water equal to soil total water
2560     
2561             do k = 1, lsoil
2562               do i = 1, im
2563                 if (slmsk(i) == 1.) then
2564                   slsoil(i,k) = smsoil(i,k)
2565                  endif
2566               enddo
2567             enddo
2568     
2569           endif   ! end if_lsm
2570     
2571     !  --- ...  total runoff is composed of drainage into water table and
2572     !           runoff at the surface and is accumulated in unit of meters
2573     
2574           if (lssav) then
2575             tem = dtf * 0.001
2576             do i = 1, im
2577               runoff(i)  = runoff(i)  + (drain(i)+runof(i)) * tem
2578               srunoff(i) = srunoff(i) + runof(i) * tem
2579     !hchaung 11/19/2007 [+1L] ;  top soil moisture unit is in mm
2580               gsoil(i)   =  gsoil(i)  + smsoil(i,1)  * dtf
2581             enddo
2582           endif
2583            
2584     !  --- ...  xw: return updated ice thickness & concentration to global array
2585     
2586           do i = 1, im
2587             if (slmsk(i) == 2.) then
2588               hice(i)  = zice(i)
2589               fice(i)  = cice(i)
2590               tisfc(i) = tice(i)
2591             else
2592               hice(i)  = 0.0
2593               fice(i)  = 0.0
2594               tisfc(i) = tsea(i)
2595             endif
2596           enddo
2597     
2598     !  --- ...  return updated smsoil and stsoil to global arrays
2599     
2600           do k = 1, lsoil
2601             do i = 1, im
2602               smc(i,k) = smsoil(i,k)
2603               stc(i,k) = stsoil(i,k)
2604               slc(i,k) = slsoil(i,k)
2605             enddo
2606           enddo
2607     
2608     !  --- ...  calculate column precipitable water "pwat"
2609     
2610           do i = 1, im
2611             pwat(i)  = 0.
2612             rqtk(i)  = 0.
2613             work2(i) = 1.0 / pgr(i)
2614           enddo
2615     
2616           do k = 1, levs
2617             do i = 1, im
2618               work1(i) = 0.0
2619             enddo
2620     
2621             if (ncld > 0) then
2622               do ic = ntcw, ntcw+ncld-1
2623                 do i = 1, im
2624                   work1(i) = work1(i) + gq0(i,k,ic)
2625                 enddo
2626               enddo
2627             endif
2628     
2629             do i = 1, im
2630               pwat(i) = pwat(i) + del(i,k)*(gq0(i,k,1)+work1(i))
2631               rqtk(i) = rqtk(i) + del(i,k)*(gq0(i,k,1)-qgrs(i,k,1))
2632             enddo
2633           enddo
2634     
2635           do i = 1, im
2636             pwat(i) = pwat(i) * (1.0/con_g)
2637             rqtk(i) = rqtk(i) * work2(i)
2638           enddo
2639     !
2640           deallocate (clw)
2641     !
2642     !     write(0,*)' pwat in gbphys',(pwat(i),i=1,im)
2643     !     if (lprnt) call mpi_quit(7)
2644     !     if (kdt == 1)  call mpi_quit(701)
2645     
2646           return
2647     !...................................
2648           end subroutine gbphys
2649     !-----------------------------------
2650     
2651