File: C:\NOAA\NEMS_11731\src\atmos\gfs\phys\gloopr.f

1            subroutine gloopr
2     !*   &    ( grid_gr,
3          &    ( grid_fld, g3d_fld,                                          &
4          &     lats_nodes_r,global_lats_r, lonsperlar, phour,               &
5          &     xlon,xlat,coszdg,COSZEN,                                     &
6          &     SLMSK,SNWDPH,SNCOVR,SNOALB,ZORL,TSEA,HPRIME,SFALB,           &
7          &     ALVSF,ALNSF ,ALVWF ,ALNWF,FACSF ,FACWF,CV,CVT ,              &
8          &     CVB  ,SWH,HLW,SFCNSW,SFCDLW,                                 &
9          &     FICE ,TISFC, SFCDSW, sfcemis,                                &
10          &     TSFLW,FLUXR, phy_f3d,slag,sdec,cdec,NBLCK,KDT,               &
11          &     global_times_r)
12     
13     !! Code Revision:
14     !! Oct 11 2009       Sarah Lu, grid_gr is replaced by grid_fld
15     !! Oct 16 2009       Sarah Lu, grid_fld%tracers used
16     !! Dec 01 2009       Sarah Lu, update fcld (instant cloud cover) in addition
17     !!                             to cldcov (cumulative cloud cover)
18     !! Dec 09 2009       Sarah Lu, (1) g3d_fld added to calling argument; (2) grrad
19     !!                   returns instant cloud cover (cldcov_v); the accumulative 
20     !!                   and instant cloud cover fields are updated after grrad call
21     !! Dec 11 2009       Sarah Lu, ldiag3d removed from grrad calling argument
22     !! Jul/Aug 2009      S. Moorthi Merged with McICA version of radiation from YuTai
23     !!
24     !!
25     !#include "f_hpm.h"
26     !
27           USE MACHINE              ,     ONLY : kind_phys, kind_grid,       &
28          &                                      kind_evod
29           USE FUNCPHYS             ,     ONLY : fpkap
30           USE PHYSCONS, fv => con_fvirt, rerth => con_rerth,
31          &              rk => con_rocp
32     
33           use module_radiation_driver,   only : radinit, grrad
34           use module_radiation_astronomy,only : astronomy
35           USE gfs_phy_tracer_config,     only : gfs_phy_tracer
36     
37           use module_radsw_parameters,  only : topfsw_type, sfcfsw_type
38           use module_radlw_parameters,  only : topflw_type, sfcflw_type
39     !
40     !! ---  for optional spectral band heating outputs
41     !!    use module_radsw_parameters,   only : NBDSW
42     !!    use module_radlw_parameters,   only : NBDLW
43     !
44           use resol_def,            ONLY: levs, levr, latr, lonr, lotgr,    &
45          &                                g_t, g_p, g_q, g_dp, g_ps,        &
46          &                                ntcw, ntoz, ncld, num_p3d,        &
47          &                                nmtvr, ntrac, levp1, nfxr, g_dpdt,&
48          &                                lgocart
49           use layout1,              ONLY: me, nodes, lats_node_r,           &
50          &                                lats_node_r_max, ipt_lats_node_r
51           use gg_def,               ONLY: coslat_r, sinlat_r
52           use date_def,             ONLY: idate
53           use namelist_physics_def, ONLY: lsswr, iaer, lslwr, sashal, ras,  &
54          &                                lssav, flgmin, ldiag3d, lggfs3d,  &
55          &                                iovr_lw, iovr_sw, isol, iems,     &
56          &                                ialb, fhlwr, fhswr, ico2, ngptc,  &
57          &                                crick_proof, norad_precip, ccnorm,&
58          &                                ictm, isubc_sw, isubc_lw, fdaer
59           use d3d_def ,             ONLY: cldcov
60           use gfs_physics_gridgr_mod, ONLY: Grid_Var_Data
61           use gfs_physics_g3d_mod,    ONLY: G3D_Var_Data
62           use mersenne_twister, only : random_setseed, random_index,        &
63          &                             random_stat
64     !
65           implicit none
66     !
67           real (kind=kind_phys), parameter :: QMIN =1.0e-10                 &
68          &,                                   Typical_pgr = 95.0            &
69          &,                                   cons0 = 0.0,  cons2 = 2.0     &
70          &,                                   pt00001=1.0e-5
71     !    &,                                   pt01=0.01
72     !
73     !  --- ...  inputs:
74           integer, intent(in) :: lats_nodes_r(nodes)
75           integer, intent(in) :: global_lats_r(latr), lonsperlar(latr)
76     
77     !*    real(kind=kind_grid) grid_gr(lonr*lats_node_r_max,lotgr)
78     
79           TYPE(Grid_Var_Data) :: grid_fld 
80           TYPE(G3D_Var_Data)  :: g3d_fld 
81     
82           integer, intent(in) :: NBLCK
83     
84     
85           real (kind=kind_phys), dimension(LONR,LATS_NODE_R), intent(in) :: &
86          &                       xlon, xlat, slmsk, snwdph, zorl, tsea,     &
87          &                       alvsf, alnsf, alvwf, alnwf, facsf, facwf,  &
88          &                       cv, cvt, cvb, FICE, tisfc, sncovr, snoalb
89     
90           real (kind=kind_phys), intent(in) ::                              &
91          &                    hprime(NMTVR,LONR,LATS_NODE_R), phour,        &
92          &                    phy_f3d(NGPTC,LEVS,NBLCK,LATS_NODE_R,NUM_P3D)
93     !
94     
95           real (kind=kind_phys), intent(inout) ::                           &
96          &                    fluxr (NFXR,LONR,LATS_NODE_R)
97     
98           integer, intent(in) :: KDT
99     !  --- ...  outputs:
100           real(kind=kind_evod), intent(out) ::                              &
101          &                    global_times_r(latr,NODES)
102     
103           real (kind=kind_phys), intent(out) ::                             &
104          &                    swh(NGPTC,LEVS,NBLCK,LATS_NODE_R),            &
105          &                    hlw(NGPTC,LEVS,NBLCK,LATS_NODE_R)
106     
107           real (kind=kind_phys),dimension(LONR,LATS_NODE_R), intent(out) :: &
108          &                    coszdg, coszen, sfcnsw, sfcdlw, tsflw,        &
109          &                    sfcdsw, SFALB, sfcemis 
110     
111           real (kind=kind_phys), intent(out) :: slag, sdec, cdec
112     
113     !! --- ...  optional spectral band heating rates
114     !!    real (kind=kind_phys), optional, intent(out) ::                   &
115     !!   &                 htrswb(NGPTC,LEVS,NBDSW,NBLCK,LATS_NODE_R),      &
116     !!   &                 htrlwb(NGPTC,LEVS,NBDLW,NBLCK,LATS_NODE_R)
117     
118     !  --- ...  locals:
119           real(kind=kind_phys) :: prsl(NGPTC,LEVS),  prslk(NGPTC,LEVS),     &
120          &                        prsi(NGPTC,LEVP1)
121     !    &                        prsi(NGPTC,LEVP1), prsik(NGPTC,LEVP1)
122     
123           real (kind=kind_phys) :: si_loc(LEVR+1)
124     
125           real (kind=kind_phys) ::                                          &
126          &                        gt(NGPTC,LEVR),                           &
127          &                        gr(NGPTC,LEVR), gr1(NGPTC,LEVR,NTRAC-1)
128     
129           real (kind=kind_phys) :: f_ice(NGPTC,LEVS), f_rain(NGPTC,LEVS),   &
130          &                         r_rime(NGPTC,LEVS)
131     
132           real (kind=kind_phys) :: cldcov_v(NGPTC,LEVS), hprime_v(NGPTC),   &
133          &                         fluxr_v(NGPTC,NFXR), vvel(NGPTC,LEVS)
134           real (kind=kind_phys) :: flgmin_v(ngptc), work1, work2
135     
136           real (kind=kind_phys) :: rinc(5), dtsw, dtlw, solcon, raddt
137     
138           real (kind=kind_phys), save :: facoz
139     
140           integer :: njeff, lon, lan, lat, iblk, lons_lat, kk
141           integer :: idat(8), jdat(8), DAYS(13), iday, imon, midmon, id
142     !     integer :: jlonr, ilan
143     
144     !  ---  variables of instantaneous calculated toa/sfc radiation fluxes
145           type (topfsw_type), dimension(NGPTC) :: topfsw
146           type (sfcfsw_type), dimension(NGPTC) :: sfcfsw
147     
148           type (topflw_type), dimension(NGPTC) :: topflw
149           type (sfcflw_type), dimension(NGPTC) :: sfcflw
150     
151     !  ---  variables used for random number generator (thread safe mode)
152           type (random_stat) :: stat
153           integer :: numrdm(LONR*LATR*2), ixseed(LONR,LATS_NODE_R,2)
154           integer :: ipseed, icsdlw(NGPTC), icsdsw(NGPTC)
155           integer, parameter :: ipsdlim = 1.0e8      ! upper limit for random seeds
156     
157     
158           integer, save :: icwp, k1oz, k2oz, midm, midp, ipsd0, iaerflg
159     
160     !  ---  number of days in a month
161           data DAYS / 31,28,31,30,31,30,31,31,30,31,30,31,30 /
162     
163     !  --- ...  control parameters: 
164     !           (some of the them may be moved into model namelist)
165     
166     !  ---  ICTM=yyyy#, controls time sensitive external data (e.g. CO2, solcon, aerosols, etc)
167     !     integer, parameter :: ICTM =   -2 ! same as ICTM=0, but add seasonal cycle
168     !                                       ! from climatology. no extrapolation.
169     !     integer, parameter :: ICTM =   -1 ! use user provided external data set for the
170     !                                       ! forecast time. no extrapolation.
171     !     integer, parameter :: ICTM =    0 ! use data at initial cond time, if not
172     !                                       ! available, use latest, no extrapolatio n.
173     !!    integer, parameter :: ICTM =    1 ! use data at the forecast time, if not
174     !                                       ! available, use latest and extrapolation.
175     !     integer, parameter :: ICTM =yyyy0 ! use yyyy data for the forecast time,
176     !                                       ! no further data extrapolation.
177     !     integer, parameter :: ICTM =yyyy1 ! use yyyy data for the fcst. if needed, do
178     !                                       ! extrapolation to match the fcst time.
179     
180     !  ---  ISOL controls solar constant data source
181     !!    integer, parameter :: ISOL  = 0  ! use prescribed solar constant
182     !     integer, parameter :: ISOL  = 1  ! use varying solar const with 11-yr cycle
183     
184     !  ---  ICO2 controls co2 data source for radiation
185     !     integer, parameter :: ICO2 = 0   ! prescribed global mean value (old opernl)
186     !!    integer, parameter :: ICO2 = 1   ! use obs co2 annual mean value only
187     !     integer, parameter :: ICO2 = 2   ! use obs co2 monthly data with 2-d variation
188     
189     !  ---  IALB controls surface albedo for sw radiation
190     !!    integer, parameter :: IALB = 0   ! use climatology alb, based on sfc type
191     !     integer, parameter :: IALB = 1   ! use modis derived alb (to be developed)
192     
193     !  ---  IEMS controls surface emissivity and sfc air/ground temp for lw radiation
194     !        ab: 2-digit control flags. a-for sfc temperature;  b-for emissivity
195     !!    integer, parameter :: IEMS = 00  ! same air/ground temp; fixed emis = 1.0
196     !!    integer, parameter :: IEMS = 01  ! same air/ground temp; varying veg typ based emis
197     !!    integer, parameter :: IEMS = 10  ! diff air/ground temp; fixed emis = 1.0
198     !!    integer, parameter :: IEMS = 11  ! diff air/ground temp; varying veg typ based emis
199     !  ---  IAER  controls aerosols scheme selections
200     !     Old definition
201     !     integer, parameter :: IAER  = 1  ! opac climatology, without volc forcing
202     !     integer, parameter :: IAER  =11  ! opac climatology, with volcanic forcing
203     !     integer, parameter :: IAER  = 2  ! gocart prognostic, without volc forcing
204     !     integer, parameter :: IAER  =12  ! gocart prognostic, with volcanic forcing
205     !     New definition in this code IAER = abc (a:volcanic; b:lw; c:sw)
206     !                             b, c values: (0:none; 1:opac; 2:gocart)
207     !  IAER =   0 --> no aerosol effect at all (volc, sw, lw)
208     !       =   1 --> only tropospheric sw aerosols, no trop-lw and volc
209     !       =  10 --> only tropospheric lw aerosols, no trop-sw and volc
210     !       =  11 --> both trop-sw and trop-lw aerosols, no volc
211     !       = 100 --> only strato-volc aeros, no trop-sw and trop-lw
212     !       = 101 --> only sw aeros (trop + volc), no lw aeros
213     !       = 110 --> only lw aeros (trop + volc), no sw aeros
214     !       = 111 --> both sw and lw aeros (trop + volc)
215     !
216     
217     !  ---  IOVR controls cloud overlapping method in radiation:
218     !     integer, parameter :: IOVR_SW = 0  ! sw: random overlap clouds
219     !!    integer, parameter :: IOVR_SW = 1  ! sw: max-random overlap clouds
220     
221     !     integer, parameter :: IOVR_LW = 0  ! lw: random overlap clouds
222     !!    integer, parameter :: IOVR_LW = 1  ! lw: max-random overlap clouds
223     
224     !  ---  ISUBC controls sub-column cloud approximation in radiation:
225     !     integer, parameter :: ISUBC_SW = 0 ! sw: without sub-col clds approx
226     !     integer, parameter :: ISUBC_SW = 1 ! sw: sub-col clds with prescribed seeds
227     !     integer, parameter :: ISUBC_SW = 2 ! sw: sub-col clds with random seeds
228     
229     !     integer, parameter :: ISUBC_LW = 0 ! lw: without sub-col clds approx
230     !     integer, parameter :: ISUBC_LW = 1 ! lw: sub-col clds with prescribed seeds
231     !     integer, parameter :: ISUBC_LW = 2 ! lw: sub-col clds with random seeds
232     
233     !  ---  iflip indicates model vertical index direction:
234     !     integer, parameter :: IFLIP = 0    ! virtical profile index from top to bottom
235           integer, parameter :: IFLIP = 1    ! virtical profile index from bottom to top
236     !
237     !    The following parameters are from gbphys
238     !
239           real (kind=kind_phys), parameter :: dxmax=-16.118095651,          &
240          &                                    dxmin=-9.800790154,           &
241          &                                    dxinv=1.0/(dxmax-dxmin)
242     
243           integer :: ierr, dimg
244           integer :: i, j, k, n, item
245     
246           logical :: change
247           logical, save :: first, sas_shal
248           data  first / .true. /
249     !
250     !  ---  for debug test use
251           real (kind=kind_phys) :: temlon, temlat, alon, alat
252           integer :: ipt
253           logical :: lprnt
254     
255     !  ---  timers:
256           real*8 :: rtc, timer1, timer2
257     !
258     !===> *** ...  begin here
259     !
260     !!
261           integer              kap,kar,kat,kau,kav,kdrlam
262           integer              ksd,ksplam,kspphi,ksq,ksr,kst
263           integer              ksu,ksv,ksz,node
264     !!
265     !     print *,' enter gloopr '
266     !
267           idat = 0
268           idat(1) = idate(4)
269           idat(2) = idate(2)
270           idat(3) = idate(3)
271           idat(5) = idate(1)
272           rinc = 0.
273     ! test repro
274     !     rinc(2) = fhour
275           rinc(2) = phour
276     !     print *,' idate ',idate
277     !     print *,' rinc ',rinc
278           call w3movdat(rinc, idat, jdat)
279     !     print *,' jdat ',jdat
280     !
281           if (ntoz .le. 0) then                ! Climatological Ozone!
282     !
283           if(me .eq. 0) WRITE (6,989) jdat(1),jdat(2),jdat(3),jdat(5)
284       989 FORMAT(' UPDATING OZONE FOR ', I4,I3,I3,I3)
285     !
286             IDAY   = jdat(3)
287             IMON   = jdat(2)
288             MIDMON = DAYS(IMON)/2 + 1
289             CHANGE = FIRST .OR.
290          &          ( (IDAY .EQ. MIDMON) .AND. (jdat(5).EQ.0) )
291     !
292             IF (CHANGE) THEN
293                 IF (IDAY .LT. MIDMON) THEN
294                    K1OZ = MOD(IMON+10,12) + 1
295                    MIDM = DAYS(K1OZ)/2 + 1
296                    K2OZ = IMON
297                    MIDP = DAYS(K1OZ) + MIDMON
298                 ELSE
299                    K1OZ = IMON
300                    MIDM = MIDMON
301                    K2OZ = MOD(IMON,12) + 1
302                    MIDP = DAYS(K2OZ)/2 + 1 + DAYS(K1OZ)
303                 ENDIF
304             ENDIF
305     !
306             IF (IDAY .LT. MIDMON) THEN
307                ID = IDAY + DAYS(K1OZ)
308             ELSE
309                ID = IDAY
310             ENDIF
311             FACOZ = real (ID-MIDM) / real (MIDP-MIDM)
312           endif
313     !
314           if (first) then
315             sas_shal = sashal .and. (.not. ras)
316     !
317             si_loc(1)=1.0
318             do k=1,levr-1
319     !*        si_loc(k+1)=si_loc(k)-grid_gr(1,g_dp+k-1)/grid_gr(1,g_ps)
320               si_loc(k+1)=si_loc(k)-grid_fld%dp(1,1,k)/grid_fld%ps(1,1) 
321             enddo
322             si_loc(levr+1)=0.0
323     
324     !  --- determin prognostic/diagnostic cloud scheme
325     
326             icwp   = 0
327             if (NTCW > 0) icwp = 1
328     
329     !  ---  generate initial permutation seed for random number generator
330     
331             if ( ISUBC_LW==2 .or. ISUBC_SW==2 ) then
332               ipsd0 = 17*idate(1) + 43*idate(2) + 37*idate(3) + 23*idate(4)
333               if ( me == 0 ) then
334                 print *,'  Radiation sub-cloud initial seed =',ipsd0,       &
335          &              ' idate =',idate
336               endif
337             endif
338     
339             iaerflg = max( mod(IAER,10), mod(IAER/10,10) ) ! flag for trop-aer scheme selection
340     
341             first = .false.
342                
343           endif         ! end_if_first
344     !
345     !===> *** ...  radiation initialization
346     !
347           dtsw  = 3600.0 * fhswr
348           dtlw  = 3600.0 * fhlwr
349     
350           raddt = min(dtsw, dtlw)
351     
352           call radinit                                                      &
353     !  ---  input:
354          &     ( si_loc, LEVR, IFLIP, idat, jdat, ICTM, ISOL, ICO2,         &
355          &       IAER, IALB, IEMS, ICWP, NUM_P3D, ISUBC_SW, ISUBC_LW,       &
356          &       IOVR_SW, IOVR_LW, me, raddt, fdaer )
357     !  ---  output: ( none )
358                                                                                                                 
359     !
360     !===> *** ...  astronomy for sw radiation calculation.
361     !
362     !     print *,' calling astronomy'
363           call astronomy                                                    &
364     !  ---  inputs:
365          &     ( lonsperlar, global_lats_r, sinlat_r, coslat_r, xlon,       &
366     !    &       fhswr, jdat, deltim,                                       &
367          &       fhswr, jdat,                                               &
368          &       LONR, LATS_NODE_R, LATR, IPT_LATS_NODE_R, lsswr, me,       &
369     !  ---  outputs:
370          &       solcon, slag, sdec, cdec, coszen, coszdg                   &
371          &      )
372     !     print *,' returned from astro'
373     
374     !
375     !===> *** ...  generate 2-d random seeds array for sub-grid cloud-radiation
376     !
377           if ( ISUBC_LW==2 .or. ISUBC_SW==2 ) then
378             ipseed = mod(nint(100.0*sqrt(phour*3600)), ipsdlim) + 1 + ipsd0
379     
380             call random_setseed                                             &
381     !  ---  inputs:
382          &     ( ipseed,                                                    &
383     !  ---  outputs:
384          &       stat                                                       &
385          &      )
386             call random_index                                               &
387     !  ---  inputs:
388          &     ( ipsdlim,                                                   &
389     !  ---  outputs:
390          &       numrdm, stat                                               &
391          &     )
392     
393             do k = 1, 2
394               do j = 1, lats_node_r
395                 lat = global_lats_r(ipt_lats_node_r-1+j)
396                 do i = 1, LONR
397                   ixseed(i,j,k) = numrdm(i+(lat-1)*LONR+(k-1)*LATR)
398                 enddo
399               enddo
400             enddo
401           endif
402     
403     
404     !
405     !===> *** ...  spectrum to grid transformation for radiation calculation.
406     !     -----------------------------------
407     !!
408     !     call f_hpmstart(61,"gr delnpe")
409     !     call f_hpmstop(61)
410     !     call f_hpmstart(62,"gr delnpo")
411     !     call f_hpmstop(62)
412     !!
413     !     call f_hpmstart(63,"gr dezouv dozeuv")
414     !
415     !     call f_hpmstop(63)
416     !!
417     !     CALL countperf(0,5,0.)
418     !     CALL synctime()
419     !     CALL countperf(1,5,0.)
420     !!
421     !     CALL countperf(0,1,0.)
422     !!
423     !     call f_hpmstart(67,"gr sumfln")
424     !     call f_hpmstop(67)
425     !!
426     !     CALL countperf(1,1,0.)
427     !     CALL countperf(0,1,0.)                                            ! hmhj
428     !
429     !     call f_hpmstart(68,"gr sumder2")                                  ! hmhj
430     !     call f_hpmstop(68)                                                ! hmhj
431     !
432           CALL countperf(1,1,0.)                                            ! hmhj
433     !
434     !
435     !===> *** ...  starting latitude loop
436     !
437           do lan=1,lats_node_r
438     !!
439              lat = global_lats_r(ipt_lats_node_r-1+lan)
440     !!
441              lons_lat = lonsperlar(lat)
442     
443     !        jlonr = (lan-1) * lonr
444     
445     !     write(0,*)' in gloopr lan=',lan,' lons_lat=',lons_lat,' lat=',lat
446     !     write (0,*)' grid_fldps=',grid_fld%ps(1:lons_lat:ngptc,lan)
447     !!
448     !$omp parallel do schedule(dynamic,1) private(lon,j,k,item,njeff,iblk,n)
449     !$omp+private(vvel,gt,gr,gr1,work1,work2,flgmin_v)
450     !$omp+private(cldcov_v,hprime_v,fluxr_v,f_ice,f_rain,r_rime)
451     !$omp+private(prslk,prsl,prsi,topfsw,sfcfsw,topflw,sfcflw)
452     !!$omp+private(prslk,prsl,prsik,prsi,topfsw,sfcfsw,topflw,sfcflw)
453     !$omp+private(icsdsw,icsdlw)
454     !!$omp+private(lprnt,ipt)
455     !$omp+private(temlon,temlat,lprnt,ipt)
456     
457     !!!$omp+private(temlon,temlat,lprnt,ipt)
458     
459             DO lon=1,lons_lat,NGPTC
460     !!
461               NJEFF = MIN(NGPTC,lons_lat-lon+1)
462               iblk  = (lon-1)/ngptc + 1
463     !
464     !  --- ...  for debug test
465     !         alon = 236.25
466     !         alat = 56.189
467               alon = 22.5
468               alat = -12.381
469               ipt = 0
470               do i = 1, njeff
471                 item = lon + i - 1
472                 temlon = xlon(item,lan) * 57.29578
473                 if (temlon < 0.0) temlon = temlon + 360.0
474                 temlat = xlat(item,lan) * 57.29578
475                 lprnt = abs(temlon-alon) < 0.5 .and. abs(temlat-alat) < 0.5
476          &          .and. kdt > 0
477                 if ( lprnt ) then
478                   ipt = i
479                   print *,' ipt=',ipt,' lon=',lon,' lan=',lan
480                   exit
481                 endif
482               enddo
483     !         lprnt = .false.
484     !!
485     !
486     
487               do i = 1, njeff
488     !           ilan = jlonr + lon+i-1
489     !           prsi(i,1) = grid_gr(ilan,g_ps)
490                 prsi(i,1) = grid_fld%ps(lon+i-1,lan)
491               enddo
492               do k = 1, LEVS
493                 do i = 1, njeff
494     !             ilan = jlonr + lon+i-1
495     !             gt(i,k)    = grid_gr(ilan,g_t+k-1)
496     !             gr(i,k)    = max(qmin,grid_gr(ilan,g_q+k-1))
497     !             prsl(i,k)  = grid_gr(ilan,g_p+k-1)
498     !             vvel(i,k)  = grid_gr(ilan,g_dpdt+k-1)
499     !             prsi(i,k+1)= prsi(i,k)-                                      &
500     !     &                    grid_gr(ilan,g_dp+k-1)
501     !*            gr(i,k)    = max(qmin,grid_fld%q(lon+i-1,lan,k))                 
502     
503                   item      = lon+i-1
504                   gt(i,k)   = grid_fld%t(item,lan,k)                           
505                   gr(i,k)   = max(qmin,grid_fld%tracers(1)%flds(item,lan,k))                 
506                   prsl(i,k)   = grid_fld%p(item,lan,k)                         
507                   vvel(i,k)   = grid_fld%dpdt(item,lan,k)                        
508                   prsi(i,k+1) = prsi(i,k) - grid_fld%dp(item,lan,k)          
509                 enddo
510               enddo
511               do i = 1, njeff
512                 prsi (i,levs+1) = 0.0
513     !           prsik(i,levs+1) = 0.0
514               enddo
515     !         do k = 1, levs
516     !           do i = 1, njeff
517     !             prslk(i,k) = (prsl(i,k)*pt00001)**rk
518     !             prsik(i,k) = (prsi(i,k)*pt00001)**rk
519     !           enddo
520     !         enddo
521     !       write (0,*)' prslk=',prslk(1,1:6),' lan=',lan
522     !
523     !       Remaining tracers
524     !
525              do n = 1, NTRAC-1
526     !           kk = g_q + n*levs
527                 do k = 1, LEVS
528                   do i = 1, njeff
529     !               ilan = jlonr + lon+i-1
530     !               gr1(i,k,n) = grid_gr(ilan,kk+k-1)
531     !*              gr1(i,k,ntoz-1) = grid_fld%oz(lon+i-1,lan,k)  
532     !*              gr1(i,k,ntcw-1) = grid_fld%cld(lon+i-1,lan,k) 
533                     gr1(i,k,n) = grid_fld%tracers(n+1)%flds(lon+i-1,lan,k)
534                   enddo
535                 enddo
536               enddo
537     
538     !!
539     !.....
540               if (levr < levs) then
541                 do i=1,njeff
542                   prsi(i,levr+1)  = prsi(i,levp1)
543                   prsl(i,levr)    = (prsi(i,levp1)+prsi(i,levr)) * 0.5
544     !             prsik(i,levr+1) = prslk(i,levp1)
545     !             prslk(i,levr)   = fpkap(prsl(i,levr))
546     !             prslk(i,levr)   = fpkap(prsl(i,levr)*1000.0)
547                 enddo
548               endif
549     !
550               if (ntoz <= 0 .or. iaerflg == 2) then
551                 do k = 1, levs
552                   do i = 1, njeff
553                     prslk(i,k) = (prsl(i,k)*pt00001)**rk
554                   enddo
555                 enddo
556               endif
557     !
558               do i=1,njeff
559                 hprime_v(i) = hprime(1,lon+i-1,lan)
560               enddo
561     !
562               do k=1,nfxr
563                 do i=1,njeff
564                   fluxr_v(i,k) = fluxr(k,lon+i-1,lan)
565                 enddo
566               enddo
567               if (NUM_P3D == 3) then
568                 do k = 1, LEVR
569                   do i = 1, njeff
570                     f_ice (i,k) = phy_f3d(i,k,iblk,lan,1)
571                     f_rain(i,k) = phy_f3d(i,k,iblk,lan,2)
572                     r_rime(i,k) = phy_f3d(i,k,iblk,lan,3)
573                   enddo
574                 enddo
575     
576                 work1 = (log(coslat_r(lat)/(lons_lat*latr)) - dxmin) * dxinv
577                 work1 = max(0.0, min(1.0,work1))
578                 work2 = flgmin(1)*work1 + flgmin(2)*(1.0-work1)
579                 do i=1,njeff
580                   flgmin_v(i) = work2
581                 enddo
582               else
583                 do i=1,njeff
584                   flgmin_v(i) = 0.0
585                 enddo
586               endif
587      
588     !  *** ...  calling radiation driver
589      
590     !
591     !     lprnt = me .eq. 0 .and. kdt .ge. 120
592     !     lprnt = me .eq. 0 .and. lan == 13
593     !     lprnt = me .eq. 0
594     !     ipt = min(91,njeff)
595           if (lprnt) then
596     !     if (kdt .gt. 85) then
597           write(0,*)' calling grrad for me=',me,' lan=',lan,' lat=',lat
598          &,' num_p3d=',num_p3d,' snoalb=',snoalb(lon,lan),' lon=',lon
599          &,' tsea=',tsea(lon,lan),' sncovr=',sncovr(lon,lan),
600          &' snwdph=',snwdph(lon,lan),' ipt=',ipt,' njeff=',njeff
601           write(0,*) ' prsi in gloopr=',prsi(ipt,1:5)
602           write(0,*) ' gt in gloopr=',gt(ipt,1:5)
603           write(0,*) ' gr in gloopr=',gr(ipt,1:5)
604           endif
605     !
606     
607     
608               call grrad
609     !  ---  inputs:
610          &     ( prsi,prsl,prslk,gt,gr,gr1,vvel,slmsk(lon,lan),             &
611          &       xlon(lon,lan),xlat(lon,lan),tsea(lon,lan),                 &
612          &       snwdph(lon,lan),sncovr(lon,lan),snoalb(lon,lan),           &
613          &       zorl(lon,lan),hprime_v,                                    &
614          &       alvsf(lon,lan),alnsf(lon,lan),alvwf(lon,lan),              &
615          &       alnwf(lon,lan),facsf(lon,lan),facwf(lon,lan),              &
616          &       fice(lon,lan),tisfc(lon,lan),                              &
617          &       solcon,coszen(lon,lan),coszdg(lon,lan),k1oz,k2oz,facoz,    &
618          &       cv(lon,lan),cvt(lon,lan),cvb(lon,lan),                     &
619          &       IOVR_SW,IOVR_LW,f_ice,f_rain,r_rime,flgmin_v,              &
620          &       icsdsw,icsdlw,NUM_P3D,NTCW-1,NCLD,NTOZ-1,NTRAC-1,NFXR,     &
621          &       dtlw,dtsw, lsswr,lslwr,lssav,sashal,norad_precip,          &
622          &       crick_proof, ccnorm,                                       &
623          &       NGPTC,njeff,LEVR,IFLIP, me, lprnt, ipt, kdt,               &
624     !  ---  outputs:
625          &       swh(1,1,iblk,lan),topfsw,sfcfsw,sfalb(lon,lan),            &
626          &       hlw(1,1,iblk,lan),topflw,sfcflw,tsflw(lon,lan),            &
627          &       sfcemis(lon,lan),cldcov_v,                                 &
628     !  ---  input/output:
629          &       fluxr_v                                                    & 
630     !! ---  optional outputs:
631     !!   &,      HTRSWB=htrswb(1,1,1,iblk,lan),                             &
632     !!   &,      HTRLWB=htrlwb(1,1,1,iblk,lan)                              &
633          &     )
634     
635     !  --- ...  radiation fluxes for other physics process or diagnostics
636     
637               if (lsswr) then
638                 do i = 1, njeff
639                   j = lon + i - 1
640                   sfcdsw(j,lan) = sfcfsw(i)%dnfxc
641                   sfcnsw(j,lan) = sfcfsw(i)%dnfxc - sfcfsw(i)%upfxc
642                 enddo
643               endif
644               if (lslwr) then
645                 do i = 1, njeff
646                   j = lon + i - 1
647                   sfcdlw(j,lan) = sfcflw(i)%dnfxc
648                 enddo
649               endif
650     
651           if (lprnt) then
652           write(0,*)' hlw=',hlw(ipt,1:10,iblk,lan),' lan=',lan,' ipt=',ipt,
653          &' njeff=',njeff
654           write(0,*)' swh=',swh(ipt,1:10,iblk,lan),' lan=',lan
655           endif
656     !
657     ! grrad routine computes cldcov_v (instant 3D cloud cover)    -- Sarah Lu
658     ! if ldiag3d is T, update cldcov (accumulative 3D cloud cover)
659     ! if lgocart is T, update fcld (instant 3D cloud cover)
660     
661               if (ldiag3d .or. lggfs3d) then
662                 do k=1,levr
663                   do i=1,njeff
664                     item = lon+i-1
665                     cldcov(k,item,lan) = cldcov(k,item,lan)                 &
666          &                                + cldcov_v(i,k) * raddt
667                   enddo
668                 enddo
669               endif
670               if (lgocart) then
671                 do k=1,levr
672                   do i=1,njeff
673                     g3d_fld%fcld(lon+i-1,lan,k) = cldcov_v(i,k) 
674                   enddo
675                 enddo
676               endif
677     
678               do k=1,nfxr
679                 do i=1,njeff
680                   fluxr(k,lon+i-1,lan) = fluxr_v(i,k)
681                 enddo
682               enddo
683               if (levr < levs) then
684                 do k=levr+1,levs
685                   do i=1,njeff
686                     hlw(i,k,iblk,lan) = hlw(i,levr,iblk,lan)
687                     swh(i,k,iblk,lan) = swh(i,levr,iblk,lan)
688                   enddo
689                 enddo
690               endif
691      
692     c$$$          write(2900+lat,*) ' ilon = ',lon
693     c$$$          write(2900+lat,'("swh",T16,"hlw")')
694     c$$$      do k=1,levs
695     c$$$         write(2900+lat,
696     c$$$     .         '(e10.3,T16,e10.3,T31,e10.3)')
697     c$$$     . swh(1,k,iblk,lan),hlw(1,k,iblk,lan)
698     c$$$       enddo
699      
700     !!
701               CALL countperf(1,12,0.)
702               ENDDO
703     !
704           enddo
705     !!
706           call f_hpmstop(69)
707     !!
708           CALL countperf(0,5,0.)
709           CALL synctime()
710           CALL countperf(1,5,0.)
711     !sela print*,'completed gloopr_v kdt=',kdt
712     !     print *,' end of gloopr '
713     !!
714           return
715           end subroutine gloopr
716