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

1     
2           SUBROUTINE SFC_NST                                                &
3     !...................................
4     !  ---  inputs:
5          &     ( IM, KM, PS, U1, V1, T1, Q1, tref, cm, ch,                  &
6          &       prsl1, prslki, slimsk, xlon, sinlat, stress,               &
7          &       sfcemis, dlwflx, sfcnsw, rain, timestep, kdt,              &
8          &       ddvel, flag_iter, flag_guess, nst_fcst, lprnt, ipr,        &
9     !  --- Input/output
10          &       tskin, tsurf, xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, &
11          &       z_c,   c_0,   c_d,   w_0, w_d, d_conv, ifd, Qrain,         &
12     !  ---  outputs:
13          &       qsurf, gflux, CMM, CHH, EVAP, HFLX, EP                     &
14          &      )
15     !
16     ! ===================================================================== !
17     !  description:                                                         !
18     !                                                                       !
19     !                                                                       !
20     !  usage:                                                               !
21     !                                                                       !
22     !    call sfc_nst                                                    !
23     !       inputs:                                                         !
24     !          ( im, km, ps, u1, v1, t1, q1, tref, cm, ch,                  !
25     !            prsl1, prslki, slimsk, xlon, sinlat, stress,               !
26     !            sfcemis, dlwflx, sfcnsw, rain, timestep, kdt,              !
27     !            ddvel, flag_iter, flag_guess, nst_fcst, lprnt, ipr,        !
28     !       input/outputs:                                                  !
29     !            tskin, tsurf, xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, !
30     !            z_c, c_0,   c_d,   w_0, w_d, d_conv, ifd, Qrain,           !
31     !  --   outputs:
32     !            qsurf, gflux, CMM, CHH, EVAP, HFLX, EP                     !
33     !           )
34     !                                                                       !
35     !                                                                       !
36     !  subprogram/functions called: w3movdat, iw3jdn, fpvs, density,        !
37     !       rhocoef, cool_skin, warm_layer, jacobi_temp.                    !
38     !                                                                       !
39     !  program history log:                                                 !
40     !         2007  -- xu li       createad original code                   !
41     !         2008  -- s. Moorthi  Adapted to the parallel version          !
42     !    may  2009  -- y.-t. hou   modified to include input lw surface     !
43     !                     emissivity from radiation. also replaced the      !
44     !                     often comfusing combined sw and lw suface         !
45     !                     flux with separate sfc net sw flux (defined       !
46     !                     as dn-up) and lw flux. added a program doc block. !
47     !    sep  2009 --  s. moorthi removed rcl and additional reformatting   !
48     !                     and optimization + made pa as input pressure unit.!
49     !         2009  -- xu li       recreatead the code                      !
50     !    feb  2010  -- s. moorthi added ass the changed made to the previous!
51     !                  version                                              !
52     !                                                                       !
53     !                                                                       !
54     !  ====================  defination of variables  ====================  !
55     !                                                                       !
56     !  inputs:                                                       size   !
57     !     im       - integer, horiz dimension                          1    !
58     !     km       - integer, vertical dimension                       1    !
59     !     ps       - real, surface pressure (Pa)                       im   !
60     !     u1, v1   - real, u/v component of surface layer wind (m/s)   im   !
61     !     t1       - real, surface layer mean temperature ( k )        im   !
62     !     q1       - real, surface layer mean specific humidity        im   !
63     !     tref     - real, reference/foundation temperature ( k )      im   !
64     !     cm       - real, surface exchange coeff for momentum (m/s)   im   !
65     !     ch       - real, surface exchange coeff heat & moisture(m/s) im   !
66     !     prsl1    - real, surface layer mean pressure (Pa)            im   !
67     !     prslki   - real,                                             im   !
68     !     slimsk   - real, sea/land/ice mask (=0/1/2)                  im   !
69     !     xlon     - real, longitude         (radians)                 im   !
70     !     sinlat   - real, sin of latitude                             im   !
71     !     stress   - real, wind stress       (n/m**2)                  im   !
72     !     sfcemis  - real, sfc lw emissivity (fraction)                im   !
73     !     dlwflx   - real, total sky sfc downward lw flux (w/m**2)     im   !
74     !     sfcnsw   - real, total sky sfc netsw flx into ocean (w/m**2) im   !
75     !     rain     - real, rainfall rate     (kg/m**2/s)               im   !
76     !     timestep - real, timestep interval (second)                  1    !
77     !     kdt      - integer, time step counter                        1    !
78     !     ddvel    - real, wind enhancement due to convection (m/s)    im   !
79     !     flag_iter- logical,                                          im   !
80     !     flag_guess-logical,                                          im   !
81     !     nst_fcst  -integer, flag 0 for no nst, 1 for uncoupled nst        !
82     !                          and 2 for coupled NST                   1    !
83     !     lprnt    - logical, control flag for check print out         1    !
84     !     ipr      - integer, grid index for check print out           1    !
85     !                                                                       !
86     !  input/outputs:
87     ! li added for oceanic components
88     !     tskin    - real, ocean surface skin temperature ( k )        im   !
89     !     tsurf    - real, ocean surface skin temperature ( k )??      im   !
90     !     xt       - real, heat content in DTL                         im   !
91     !     xs       - real, salinity  content in DTL                    im   !
92     !     xu       - real, u-current content in DTL                    im   !
93     !     xv       - real, v-current content in DTL                    im   !
94     !     xz       - real, DTL thickness                               im   !
95     !     zm       - real, MXL thickness                               im   !
96     !     xtts     - real, d(xt)/d(ts)                                 im   !
97     !     xzts     - real, d(xz)/d(ts)                                 im   !
98     !     dt_cool  - real, Sub-layer cooling amount                    im   !
99     !     d_conv   - real, thickness of Free Convection Layer (FCL)    im   !
100     !     z_c      - Sub-layer cooling thickness                       im   !
101     !     c_0      - coefficient1 to calculate d(Tz)/d(Ts)             im   !
102     !     c_d      - coefficient2 to calculate d(Tz)/d(Ts)             im   !
103     !     w_0      - coefficient3 to calculate d(Tz)/d(Ts)             im   !
104     !     w_d      - coefficient4 to calculate d(Tz)/d(Ts)             im   !
105     !     ifd      - real, index to start DTM run or not               im   !
106     !     Qrain    - real, sensible heat flux due to rainfall (watts)  im   !
107     
108     !  outputs:                                                             !
109     
110     !     qsurf    - real, surface air saturation specific humidity    im   !
111     !     gflux    - real, soil heat flux (w/m**2)                     im   !
112     !     cmm      - real,                                             im   !
113     !     chh      - real,                                             im   !
114     !     evap     - real, evaperation from latent heat flux           im   !
115     !     hflx     - real, sensible heat flux                          im   !
116     !     ep       - real, potential evaporation                       im   !
117     !                                                                       !
118     ! ===================================================================== !
119           USE MACHINE , ONLY : kind_phys
120           USE FUNCPHYS, ONLY : fpvs
121           USE PHYSCONS, HVAP => con_HVAP
122          &,             CP => con_CP, HFUS => con_HFUS, JCAL => con_JCAL
123          &,             EPS => con_eps, EPSM1 => con_epsm1
124          &,             RVRDM1 => con_FVirt, RD => con_RD
125          &,             RHW0 => con_rhw0,SBC => con_sbc,pi => con_pi
126           USE date_def, only: idate
127           USE module_nst_parameters, ONLY : t0K,cp_w,omg_m,omg_sh,
128          &    sigma_r,gray,solar_time_6am,Ri_c,z_w_max,delz,wd_max,
129          &    rad2deg,const_rot,tau_min,tw_max,sst_max
130           USE module_nst_water_prop, ONLY: solar_time_from_julian,
131          &     density,rhocoef,compjd,grv 
132          &,    sw_ps_9b
133           USE nst_module, ONLY : cool_skin,dtm_1p,cal_w,cal_ttop,convdepth,
134          &  dtm_1p_fca,dtm_1p_tla,dtm_1p_mwa,dtm_1p_mda,dtm_1p_mta,dtl_reset
135     !
136           implicit none
137     !
138     !  ---  constant parameters:
139           real (kind=kind_phys), parameter :: cpinv=1.0/cp, HVAPI=1.0/HVAP
140           real (kind=kind_phys), parameter :: EMISSIV=0.97
141           real (kind=kind_phys), parameter :: f24   = 24.0     ! hours/day
142           real (kind=kind_phys), parameter :: f1440 = 1440.0   ! minutes/day
143     
144     !  ---  inputs:
145           integer, intent(in) :: im, km, kdt, ipr, nst_fcst
146           real (kind=kind_phys), dimension(im), intent(in) :: ps, u1, v1,   &
147          &       t1, q1, tref, cm, ch, prsl1, prslki, slimsk, xlon,         &
148          &       sinlat, stress, sfcemis, dlwflx, sfcnsw, rain, ddvel
149     
150           real (kind=kind_phys), intent(in) :: timestep
151     
152           logical, intent(in) :: flag_iter(im), flag_guess(im), lprnt
153     
154     !  ---  input/outputs:
155     ! control variables of DTL system (5+2) and SL (2) and coefficients for d(Tz)/d(Ts) calculation
156           real (kind=kind_phys), dimension(im), intent(inout) :: tskin,     &
157          &      tsurf, xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool,         &
158          &      z_c, c_0, c_d, w_0, w_d, d_conv, ifd, qrain
159     
160     !  ---  outputs:
161           real (kind=kind_phys), dimension(im), intent(out) ::              &
162          &       qsurf, gflux, cmm, chh, evap, hflx, ep
163     
164     !
165     !     Locals
166     !
167           integer :: k,i,iter
168     !
169           real (kind=kind_phys), dimension(im) ::  Q0, QSS, RCH, 
170          &                     rho_a, THETA1, TV1, WIND, wndmag
171     
172           real(kind=kind_phys) elocp,tem
173     !
174     !    NSTM related prognostic fields
175     !
176           logical flag(im)
177           real (kind=kind_phys), dimension(im) ::        
178          &   xt_old, xs_old, xu_old, xv_old, xz_old,zm_old,xtts_old,
179          &   xzts_old, ifd_old, tref_old, tskin_old, dt_cool_old,z_c_old
180     
181           real(kind=kind_phys) ulwflx(im), nswsfc(im)
182     !     real(kind=kind_phys) rig(im),
183     !    &                     ulwflx(im),dlwflx(im),
184     !    &                     slrad(im),nswsfc(im)
185           real(kind=kind_phys) alpha,beta,rho_w,F_nsol,sss,sep,soltim,
186          & cosa,sina,taux,tauy,grav,dz,t0,ttop0,ttop
187     
188           integer :: iyear,imon,iday,ihr,imin,jd
189           integer :: idat(8),jdat(8)
190     
191           real(kind=kind_phys) Le,fc,dwat,dtmp,wetc,alfac,ustar_a,rich
192           real(kind=kind_phys) Rnl_Ts,Hs_Ts,Hl_Ts,RF_Ts,Q_Ts
193           real(kind=kind_phys) jday,jday_old,fw,Q_warm
194           real(kind=kind_phys) t12,alon,alat,es,Qs,tsea,sstc,dta
195           real(kind=kind_phys) fjd1,jd1,jd0
196           real(kind=kind_phys) rinc(5) 
197     
198     !  external functions called: iw3jdn
199           integer :: iw3jdn
200     !======================================================================================================
201     cc
202           PARAMETER (elocp=hvap/cp)
203     
204           sss = 34.0             ! temporarily, when sea surface salinity data is not ready
205     
206           idat(1) = idate(4)
207           idat(2) = idate(2)
208           idat(3) = idate(3)
209           idat(4) = 0
210           idat(5) = idate(1)
211           idat(6) = 0
212           idat(7) = 0
213           idat(8) = 0
214     
215           rinc(1) = 0.
216           rinc(2) = 0.
217           rinc(3) = float(kdt)*timestep/60.0
218           rinc(4) = 0.
219           rinc(5) = 0.
220           call w3movdat(rinc, idat, jdat)
221     
222           iyear = jdat(1)
223           imon  = jdat(2)
224           iday  = jdat(3)
225           ihr   = jdat(5)
226           imin  = jdat(6)
227     
228     !  --- ...  calculate forecast julian day and fraction of julian day
229     
230           jd0 = iw3jdn(1899,12,31)
231           jd1 = iw3jdn(iyear,imon,iday)
232     
233     !  --- ...  unlike in normal applications, where day starts from 0 hr,
234     !           in astronomy applications, day stats from noon.
235           if (ihr < 12) then
236             jd1 = jd1 - 1
237             fjd1= 0.5 + float(ihr)/f24 + float(imin)/f1440
238           else
239             fjd1= float(ihr - 12)/f24 + float(imin)/f1440
240           endif
241           jday  = jd1 - jd0 + fjd1 
242     !
243     ! FLAG for open water
244     !
245           DO I = 1, IM
246              FLAG(I) = SLIMSK(I).EQ. 0. .AND. flag_iter(i)
247           ENDDO
248     
249     !
250     !  save nst-related prognostic fields for guess run
251     !
252           do i=1, im
253             if((SLIMSK(I).EQ. 0.) .AND. flag_guess(i)) then
254               xt_old(i)      = xt(i)
255               xs_old(i)      = xs(i)
256               xu_old(i)      = xu(i)
257               xv_old(i)      = xv(i)
258               xz_old(i)      = xz(i)
259               zm_old(i)      = zm(i)
260               xtts_old(i)    = xtts(i)
261               xzts_old(i)    = xzts(i)
262               ifd_old(i)     = ifd(i)
263               tskin_old(i)   = tskin(i)
264               dt_cool_old(i) = dt_cool(i)
265               z_c_old(i)     = z_c(i)
266             endif
267           enddo
268     
269     !  --- ...  initialize variables. all units are m.k.s. unless specified.
270     !           ps is in pascals, wind is wind speed, theta1 is surface air
271     !           estimated from level 1 temperature, rho_a is air density and
272     !           qss is saturation specific humidity at the water surface
273     !!
274           do i = 1, im
275             if ( flag(i) ) then
276     
277               nswsfc(i) = sfcnsw(i) ! net solar radiation at the air-sea surface (positive downward)
278               wndmag(i) = sqrt(u1(i)*u1(i) + v1(i)*v1(i))
279               wind(i)   = wndmag(i) + max( 0.0, min( ddvel(i), 30.0 ) )
280               wind(i)   = max( wind(i), 1.0 )
281     
282               q0(i)     = max(q1(i), 1.0e-8)
283               theta1(i) = t1(i) * prslki(i)
284               tv1(i)    = t1(i) * (1.0 + rvrdm1*q0(i))
285               rho_a(i)  = prsl1(i) / (rd*tv1(i))
286               qss(i)    = fpvs(tsurf(i))                          ! pa
287               qss(i)    = eps*qss(i) / (ps(i) + epsm1*qss(i))     ! pa
288     !
289               evap(i)    = 0.0
290               hflx(i)    = 0.0
291               gflux(i)   = 0.0
292               ep(i)      = 0.0
293     
294     !  --- ...  rcp = rho cp ch v
295     
296               rch(i)     = rho_a(i) * cp * ch(i) * wind(i)
297               cmm(i)     = cm (i)   * wind(i)
298               chh(i)     = rho_a(i) * ch(i) * wind(i)
299     
300     !  --- ...  latent and sensible heat flux over open water with tskin
301     !           at previous time step
302               evap(i)    = elocp * rch(i) * (qss(i) - q0(i))
303               qsurf(i)   = qss(i)
304               hflx(i)    = rch(i) * (tsurf(i) - theta1(i))
305     
306     !     if (lprnt .and. i == ipr) print *,' tskin=',tskin(i),' theta1=',
307     !    & theta1(i),' hflx=',hflx(i),' t1=',t1(i),'prslki=',prslki(i)
308     !    &,' tsurf=',tsurf(i)
309             endif
310           enddo
311     
312     ! Run NST Model: DTM + SLM   
313     ! 
314           do i = 1, im
315             if ( flag(i) ) then
316               tsea      = tsurf(i)
317               t12       = tsea*tsea
318               ulwflx(i) = sfcemis(i) * sbc * t12 * t12
319               alon      = xlon(i)*rad2deg
320               grav      = grv(sinlat(i))
321               CALL solar_time_from_julian(jday,alon,soltim)
322               CALL density(tsea,SSS,rho_w)                     ! Sea water density
323               CALL rhocoef(tsea,SSS,rho_w,alpha,beta)          ! alpha & beta
324     !
325     !  CALCULATE SENSIBLE HEAT FLUX due to rainfall 
326     !
327               Le       = (2.501-.00237*tsea)*1e6
328               dwat     = 2.11e-5*(T1(I)/t0K)**1.94               ! water vapor diffusivity
329               dtmp     = (1.+3.309e-3*(T1(I)-t0K)-1.44e-6*(T1(I)-t0K)*
330          &              (T1(I)-t0K))*0.02411/(rho_a(I)*CP)       ! heat diffusivity
331               wetc     = 622.0*Le*QSS(I)/(rd*t1(i)*t1(i))
332               alfac    = 1/(1+(wetc*Le*dwat)/(CP*dtmp))          ! wet bulb factor
333               Qrain(I) =  (1000.*rain(I)/rho_w)*alfac*cp_w*
334          &                (tsea-T1(I)+(1000.*QSS(I)-1000.*Q0(I))*Le/CP)
335     
336     !  --- ...  input non solar heat flux as upward = positive to models here
337     
338               f_nsol = hflx(i) + evap(i) + ulwflx(i) - dlwflx(i)
339          &           + omg_sh*Qrain(i)
340     
341     !     if (lprnt .and. i == ipr) print *,' f_nsol=',f_nsol,' hflx=',
342     !    &hflx(i),' evap=',evap(i),' ulwflx=',ulwflx(i),' dlwflx=',dlwflx(i)
343     !    &,' omg_sh=',omg_sh,' qrain=',qrain(i)
344     
345               sep      =  sss*(evap(i)/Le-rain(i))/rho_w
346               ustar_a  = sqrt(stress(i)/rho_a(i))          ! air friction velocity
347     !
348     !  Sensitivities of heat flux components to Ts
349     !
350               Rnl_Ts = 4.0*gray*sigma_r*tsea*tsea*tsea     ! d(Rnl)/d(Ts)
351               Hs_Ts = rch(i)
352               Hl_Ts = rch(i)*elocp*eps*hvap*qss(i)/(rd*t12)
353               RF_Ts = (1000.*rain(i)/rho_w)*alfac*cp_w*(1.0+rch(i)*Hl_Ts)
354               Q_Ts  = Rnl_Ts + Hs_Ts + Hl_Ts + omg_sh*RF_Ts
355     !
356     !    run sub-layer cooling model/parameterization & calculate c_0, c_d
357     !
358               call cool_skin(ustar_a,F_nsol,NSWSFC(i),evap(i),sss,alpha,beta
359          &,   rho_w,rho_a(i),tsea,Q_Ts,Hl_Ts,grav,Le
360          &,   dt_cool(i),z_c(i),c_0(i),c_d(i))
361     
362               tem  = 1.0 / wndmag(i)
363               cosa = u1(i)*tem
364               sina = v1(i)*tem
365               taux = max(stress(i),tau_min)*cosa
366               tauy = max(stress(i),tau_min)*sina
367               fc   = const_rot*sinlat(i)
368     !
369     !     Run DTM-1p system
370     !
371               if ( (soltim > solar_time_6am .and. ifd(i) == 0.0) ) then
372               else
373                 ifd(i) = 1.0
374     !
375     !     Calculate FCL thickness with current forcing and previous time's profile
376     !
377     !     if (lprnt .and. i == ipr) print *,' beg xz=',xz(i)
378     
379                 if ( F_nsol > 0.0 .and. xt(i) > 0.0 ) then
380                   call convdepth(kdt,timestep,NSWSFC(i),F_nsol,sss,sep,rho_w
381          &,            alpha,beta,xt(i),xs(i),xz(i),d_conv(i))
382                 else
383                   d_conv(i) = 0.0
384                 endif
385     
386     !     if (lprnt .and. i == ipr) print *,' beg xz1=',xz(i)
387     !
388     !    determine rich: wind speed dependent (right now)
389     !
390     !           if ( wind(i) < 1.0 ) then
391     !             rich = 0.25 + 0.03*wind(i)
392     !           elseif ( wind(i) >= 1.0 .and. wind(i) < 1.5 ) then
393     !             rich = 0.25 + 0.1*wind(i)
394     !           elseif ( wind(i) >= 1.5 .and. wind(i) < 6.0 ) then
395     !             rich = 0.25 + 0.6*wind(i)
396     !           elseif ( wind(i) >= 6.0 ) then
397     !             rich = 0.25 + min(0.8*wind(i),0.50)
398     !           endif
399     
400                 rich = Ri_c
401              
402                 call dtm_1p(kdt,timestep,rich,taux,tauy,NSWSFC(i),
403          &      F_nsol,sss,sep,Q_Ts,Hl_Ts,rho_w,alpha,beta,alon,sinlat(i),
404          &      soltim,grav,Le,d_conv(i),
405          &      xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i))
406     
407     !     if (lprnt .and. i == ipr) print *,' beg xz2=',xz(i)
408     
409     !  Apply MDA
410                 if ( xt(i) > 0.0 ) then
411                   call dtm_1p_mda(xt(i),xtts(i),xz(i),xzts(i))
412                   if ( xz(i) >= z_w_max ) then
413                     call dtl_reset(xt(i),xs(i),xu(i),xv(i),xz(i),xtts(i),
414          &                                                       xzts(i))
415     
416     !     if (lprnt .and. i == ipr) print *,' beg xz3=',xz(i),' z_w_max='
417     !    &,z_w_max
418                   endif
419     
420     !  Apply FCA
421                   if ( d_conv(i) > 0.0 ) then
422                     call dtm_1p_fca(d_conv(i),xt(i),xtts(i),xz(i),xzts(i))
423                     if ( xz(i) >= z_w_max ) then
424                       call dtl_reset
425          &              (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i))
426                     endif
427                   endif
428     
429     !     if (lprnt .and. i == ipr) print *,' beg xz4=',xz(i)
430     
431     !  Apply TLA
432                   dz = min(xz(i),max(d_conv(i),delz))
433     !
434                   call sw_ps_9b(delz,fw)
435                   Q_warm = fw*NSWSFC(i)-F_nsol    !total heat absorbed in warm layer
436                   if ( Q_warm > 0.0 ) then
437                     call cal_ttop(kdt,timestep,Q_warm,rho_w,dz,
438          &                        xt(i),xz(i),ttop0)
439     
440     !     if (lprnt .and. i == ipr) print *,' d_conv=',d_conv(i),' delz=',
441     !    &delz,' kdt=',kdt,' timestep=',timestep,' nswsfc=',nswsfc(i),
442     !    &' f_nsol=',f_nsol,' rho_w=',rho_w,' dz=',dz,' xt=',xt(i),
443     !    &' xz=',xz(i),' qrain=',qrain(i)
444     
445                     ttop = ((xt(i)+xt(i))/xz(i))*(1.0-dz/((xz(i)+xz(i))))
446     
447     !     if (lprnt .and. i == ipr) print *,' beg xz4a=',xz(i)
448     !    &,' ttop=',ttop,' ttop0=',ttop0,' xt=',xt(i),' dz=',dz
449     !    &,' xznew=',(xt(i)+sqrt(xt(i)*(xt(i)-dz*ttop0)))/ttop0
450     
451                     if ( ttop > ttop0 ) then
452                       call dtm_1p_tla(dz,ttop0,xt(i),xtts(i),xz(i),xzts(i))
453     
454     !     if (lprnt .and. i == ipr) print *,' beg xz4b=',xz(i),'z_w_max=',
455     !    &z_w_max
456                       if ( xz(i) >= z_w_max ) then
457                         call dtl_reset
458          &                   (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i))
459                       endif
460                     endif
461                   endif           ! if ( Q_warm > 0.0 ) then
462     
463     !     if (lprnt .and. i == ipr) print *,' beg xz5=',xz(i)
464     
465     !  Apply MWA
466                   t0 = (xt(i)+xt(i))/xz(i)
467                   if ( t0 > tw_max ) then
468                     call dtm_1p_mwa(xt(i),xtts(i),xz(i),xzts(i))
469                     if ( xz(i) >= z_w_max ) then
470                       call dtl_reset
471          &                 (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i))
472                     endif
473                   endif
474     
475     !     if (lprnt .and. i == ipr) print *,' beg xz6=',xz(i)
476     
477     !  Apply MTA
478                   sstc = tref(i) + (xt(i)+xt(i))/xz(i) - dt_cool(i)
479                   if ( sstc > sst_max ) then
480                     dta = sstc - sst_max
481                     call  dtm_1p_mta(dta,xt(i),xtts(i),xz(i),xzts(i))
482     !               write(*,'(a,F3.0,7F8.3)') 'MTA, sstc,dta :',slimsk(i),
483     !    &          sstc,dta,tref(i),xt(i),xz(i),2.0*xt(i)/xz(i),dt_cool(i)
484                    if ( xz(i) >= z_w_max ) then
485                       call dtl_reset
486          &                 (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i))
487                     endif
488                   endif
489     !
490                 endif             ! if ( xt(i) > 0.0 ) then
491     
492               endif             ! IF (solar_time > solar_time_6am .and. ifd(i) == 0.0 ) THEN: too late to start the first day
493     
494     !     if (lprnt .and. i == ipr) print *,' beg xz7=',xz(i)
495     
496               tsurf(i) = max(271.0, tref(i)+(xt(i)+xt(i))/xz(i)-dt_cool(i))
497     
498     !     if (lprnt .and. i == ipr) print *,' tsurf=',tsurf(i),' tref=',
499     !    &tref(i),' xz=',xz(i),' dt_cool=',dt_cool(i)
500     
501               if ( xt(i) > 0.0 ) then
502                 call cal_w(kdt,xz(i),xt(i),xzts(i),xtts(i),w_0(i),w_d(i))
503               else
504                 w_0(i) = 0.0
505                 w_d(i) = 0.0
506               endif
507     
508     !         if ( xt(i) > 0.0 ) then
509     !           rig(i) = grav*xz(i)*xz(i)*(alpha*xt(i)-beta*xs(i))
510     !    &             /(2.0*(xu(i)*xu(i)+xv(i)*xv(i)))
511     !         else
512     !           rig(i) = 0.25
513     !         endif
514     
515     !         Qrain(i) = rig(i)
516               zm(i) = wind(i)
517             
518             ENDIF
519           ENDDO
520     
521     ! restore NST-related prognostic fields for guess run
522           do i=1, im
523             if((slimsk(I) == 0.) ) then
524               if(flag_guess(i)) then
525                 xt(i)      = xt_old(i)
526                 xs(i)      = xs_old(i)
527                 xu(i)      = xu_old(i)
528                 xv(i)      = xv_old(i)
529                 xz(i)      = xz_old(i)
530                 zm(i)      = zm_old(i)
531                 xtts(i)    = xtts_old(i)
532                 xzts(i)    = xzts_old(i)
533                 ifd(i)     = ifd_old(i)
534                 tskin(i)   = tskin_old(i)
535                 dt_cool(i) = dt_cool_old(i)
536                 z_c(i)     = z_c_old(i)
537               else
538                 if ( nst_fcst > 1 ) then
539                   tskin(i) = tsurf(i)
540                 endif                  ! if ( nst_fcst > 1 ) then
541               endif                    ! if(flag_guess(i)) then
542             endif                      ! if((SLIMSK(I).EQ. 0.) ) then
543           enddo
544     
545     !     if (lprnt .and. i == ipr) print *,' beg xz8=',xz(i)
546     
547           if ( nst_fcst > 1 ) then
548     !  --- ...  latent and sensible heat flux over open water with updated tskin
549             do i = 1, im
550               if ( flag(i) ) then
551                 qss(i)   = fpvs( tskin(i) )
552                 qss(i)   = eps*qss(i) / (ps(i) + epsm1*qss(i))
553                 qsurf(i) = qss(i)
554                 evap(i)  = elocp*rch(i) * (qss(i) - q0(i))
555                 hflx(i)  = rch(i) * (tskin(i) - theta1(i))
556               endif
557             enddo
558           endif                   ! if ( nst_fcst > 1 ) then
559     
560     !
561           do i=1,im
562             if ( flag(i) ) then
563               tem     = 1.0 / rho_a(i)
564               hflx(i) = hflx(i) * tem * cpinv
565               evap(i) = evap(i) * tem * hvapi
566             endif
567           enddo
568     !
569     !     if (lprnt) print *,' tskin=',tskin(ipr)
570     
571           return
572           end
573