/[MITgcm]/MITgcm/pkg/exf/exf_getffields.F
ViewVC logotype

Diff of /MITgcm/pkg/exf/exf_getffields.F

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

revision 1.2 by heimbach, Mon Jul 30 20:41:20 2001 UTC revision 1.3 by heimbach, Fri Jan 11 19:24:24 2002 UTC
# Line 2  c $Header$ Line 2  c $Header$
2    
3  #include "EXF_CPPOPTIONS.h"  #include "EXF_CPPOPTIONS.h"
4    
5        subroutine exf_GetFFields(        subroutine exf_GetFFields( mycurrenttime, mycurrentiter, mythid )
      I                           mycurrenttime,  
      I                           mycurrentiter,  
      I                           mythid  
      &                         )  
6    
7  c     ==================================================================  c     ==================================================================
8  c     SUBROUTINE exf_GetFFields  c     SUBROUTINE exf_GetFFields
# Line 137  c     == global variables == Line 133  c     == global variables ==
133  #include "exf_fields.h"  #include "exf_fields.h"
134  #include "exf_constants.h"  #include "exf_constants.h"
135    
136    #ifdef ALLOW_AUTODIFF_TAMC
137    #include "tamc.h"
138    #endif
139    
140  c     == routine arguments ==  c     == routine arguments ==
141    
142        integer mythid        integer mythid
# Line 165  c     == local variables == Line 165  c     == local variables ==
165        _RL     re        _RL     re
166        _RL     rdn        _RL     rdn
167        _RL     rh        _RL     rh
168          _RL     ssttmp
169        _RL     ssq        _RL     ssq
170        _RL     stable        _RL     stable
       _RL     tau  
171        _RL     tstar        _RL     tstar
172        _RL     t0        _RL     t0
173        _RL     ustar        _RL     ustar
174        _RL     uzn        _RL     uzn
175          _RL     shn
176        _RL     xsq        _RL     xsq
177        _RL     x        _RL     x
178          _RL     tau
179        _RL     evap(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL     evap(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
180    #ifdef ALLOW_AUTODIFF_TAMC
181          integer ikey_1
182          integer ikey_2
183    #endif
184  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
185    
186          _RL     ustmp
187        _RL     us        _RL     us
188        _RL     cw        _RL     cw
189        _RL     sw        _RL     sw
# Line 207  c     determine forcing field records Line 214  c     determine forcing field records
214    
215  #ifdef ALLOW_BULKFORMULAE  #ifdef ALLOW_BULKFORMULAE
216    
217    #if (defined (ALLOW_ATM_TEMP) || defined (ALLOW_ATM_WIND))
218          aln = log(ht/zref)
219    #endif
220    
221  c     Determine where we are in time and set counters, flags and  c     Determine where we are in time and set counters, flags and
222  c     the linear interpolation factors accordingly.  c     the linear interpolation factors accordingly.
223  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
224  c     Atmospheric temperature.  c     Atmospheric temperature.
225        call exf_set_atemp( atemp        call exf_set_atemp( mycurrenttime, mycurrentiter, mythid )
      &                  , mycurrenttime, mycurrentiter, mythid )  
226    
227  c     Atmospheric humidity.  c     Atmospheric humidity.
228        call exf_set_aqh( aqh        call exf_set_aqh( mycurrenttime, mycurrentiter, mythid )
      &                , mycurrenttime, mycurrentiter, mythid )  
229    
230  c     Net long wave radiative flux.  c     Net long wave radiative flux.
231        call exf_set_lwflux( lwflux        call exf_set_lwflux( mycurrenttime, mycurrentiter, mythid )
      &                   , mycurrenttime, mycurrentiter, mythid )  
232    
233  c     Net short wave radiative flux.  c     Net short wave radiative flux.
234        call exf_set_swflux( swflux        call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )
      &                   , mycurrenttime, mycurrentiter, mythid )  
235    
236  c     Precipitation.  c     Precipitation.
237        call exf_set_precip( precip        call exf_set_precip( mycurrenttime, mycurrentiter, mythid )
      &                   , mycurrenttime, mycurrentiter, mythid )  
238    
239        aln = log(ht/zref)  #ifdef ALLOW_ATEMP_CONTROL
240          call ctrl_getatemp ( mycurrenttime, mycurrentiter, mythid )
241    #endif
242    
243    #ifdef ALLOW_AQH_CONTROL
244          call ctrl_getaqh ( mycurrenttime, mycurrentiter, mythid )
245    #endif
246    
247  #else  #else
248    
249  c     Atmospheric heat flux.  c     Atmospheric heat flux.
250        call exf_set_hflux( hflux        call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )
      &                  , mycurrenttime, mycurrentiter, mythid )  
251    
252  c     Salt flux.  c     Salt flux.
253        call exf_set_sflux( sflux        call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )
      &                  , mycurrenttime, mycurrentiter, mythid )  
254    
255  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
256  c     Net short wave radiative flux.  c     Net short wave radiative flux.
257        call exf_set_swflux( swflux        call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )
      &                   , mycurrenttime, mycurrentiter, mythid )  
258    
259  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
260    
# Line 253  c     Net short wave radiative flux. Line 262  c     Net short wave radiative flux.
262    
263  #ifdef ALLOW_ATM_WIND  #ifdef ALLOW_ATM_WIND
264  c     Zonal wind.  c     Zonal wind.
265        call exf_set_uwind( uwind        call exf_set_uwind( mycurrenttime, mycurrentiter, mythid )
      &                  , mycurrenttime, mycurrentiter, mythid )  
266    
267  c     Meridional wind.  c     Meridional wind.
268        call exf_set_vwind( vwind        call exf_set_vwind( mycurrenttime, mycurrentiter, mythid )
269       &                  , mycurrenttime, mycurrentiter, mythid )  
270    #ifdef ALLOW_UWIND_CONTROL
271          call ctrl_getuwind ( mycurrenttime, mycurrentiter, mythid )
272    #endif
273    
274    #ifdef ALLOW_VWIND_CONTROL
275          call ctrl_getvwind ( mycurrenttime, mycurrentiter, mythid )
276    #endif
277    
278  #else  #else
279  c     Zonal wind stress.  c     Zonal wind stress.
280        call exf_set_ustress( ustress        call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )
      &                    , mycurrenttime, mycurrentiter, mythid )  
281    
282  c     Meridional wind stress.  c     Meridional wind stress.
283        call exf_set_vstress( vstress        call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )
      &                    , mycurrenttime, mycurrentiter, mythid )  
284    
285  #endif /* ALLOW_ATM_WIND */  #endif /* ALLOW_ATM_WIND */
286    
287  #else /* ALLOW_BULKFORMULAE undefined */  #else /* ALLOW_BULKFORMULAE undefined */
288  c     Atmospheric heat flux.  c     Atmospheric heat flux.
289        call exf_set_hflux( hflux        call exf_set_hflux( mycurrenttime, mycurrentiter, mythid )
      &                  , mycurrenttime, mycurrentiter, mythid )  
290    
291  c     Salt flux.  c     Salt flux.
292        call exf_set_sflux( sflux        call exf_set_sflux( mycurrenttime, mycurrentiter, mythid )
      &                  , mycurrenttime, mycurrentiter, mythid )  
293    
294  c     Zonal wind stress.  c     Zonal wind stress.
295        call exf_set_ustress( ustress        call exf_set_ustress( mycurrenttime, mycurrentiter, mythid )
      &                    , mycurrenttime, mycurrentiter, mythid )  
296    
297  c     Meridional wind stress.  c     Meridional wind stress.
298        call exf_set_vstress( vstress        call exf_set_vstress( mycurrenttime, mycurrentiter, mythid )
      &                    , mycurrenttime, mycurrentiter, mythid )  
299    
300  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
301  c     Net short wave radiative flux.  c     Net short wave radiative flux.
302        call exf_set_swflux( swflux        call exf_set_swflux( mycurrenttime, mycurrentiter, mythid )
303       &                   , mycurrenttime, mycurrentiter, mythid )  #endif
 #endif /* ALLOW_KPP */  
304    
305  #endif /* ALLOW_BULKFORMULAE */  #endif /* ALLOW_BULKFORMULAE */
306    
307  c     Loop over tiles.  c     Loop over tiles.
308    #ifdef ALLOW_AUTODIFF_TAMC
309    C--   HPF directive to help TAMC
310    CHPF$ INDEPENDENT
311    #endif /* ALLOW_AUTODIFF_TAMC */
312        do bj = mybylo(mythid),mybyhi(mythid)        do bj = mybylo(mythid),mybyhi(mythid)
313    #ifdef ALLOW_AUTODIFF_TAMC
314    C--    HPF directive to help TAMC
315    CHPF$  INDEPENDENT
316    #endif
317          do bi = mybxlo(mythid),mybxhi(mythid)          do bi = mybxlo(mythid),mybxhi(mythid)
318    
319            k = 1            k = 1
320    
321            do j = 1-oly,sny+oly            do j = 1-oly,sny+oly
# Line 306  c     Loop over tiles. Line 323  c     Loop over tiles.
323    
324  #ifdef ALLOW_BULKFORMULAE  #ifdef ALLOW_BULKFORMULAE
325    
326    #ifdef ALLOW_AUTODIFF_TAMC
327                   act1 = bi - myBxLo(myThid)
328                   max1 = myBxHi(myThid) - myBxLo(myThid) + 1
329                   act2 = bj - myByLo(myThid)
330                   max2 = myByHi(myThid) - myByLo(myThid) + 1
331                   act3 = myThid - 1
332                   max3 = nTx*nTy
333                   act4 = ikey_dynamics - 1
334    
335                   ikey_1 = i
336         &              + sNx*(j-1)
337         &              + sNx*sNy*act1
338         &              + sNx*sNy*max1*act2
339         &              + sNx*sNy*max1*max2*act3
340         &              + sNx*sNy*max1*max2*max3*act4
341    #endif
342    
343  c             Compute the turbulent surface fluxes.  c             Compute the turbulent surface fluxes.
344  c                  (Bulk formulae estimates)  c                  (Bulk formulae estimates)
345    
346  #ifdef ALLOW_ATM_WIND  #ifdef ALLOW_ATM_WIND
347  c             Wind speed and direction.  c             Wind speed and direction.
348                us = sqrt(uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +                ustmp = uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +
349       &                  vwind(i,j,bi,bj)*vwind(i,j,bi,bj))       &                vwind(i,j,bi,bj)*vwind(i,j,bi,bj)
350                if ( us .ne. 0. _d 0 ) then                if ( ustmp .ne. 0. _d 0 ) then
351                    us = sqrt(ustmp)
352                  cw = uwind(i,j,bi,bj)/us                  cw = uwind(i,j,bi,bj)/us
353                  sw = vwind(i,j,bi,bj)/us                  sw = vwind(i,j,bi,bj)/us
354                else                else
355                    us = 0. _d 0
356                  cw = 0. _d 0                  cw = 0. _d 0
357                  sw = 0. _d 0                  sw = 0. _d 0
358                endif                endif
# Line 324  c             Wind speed and direction. Line 360  c             Wind speed and direction.
360  #else  #else
361  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
362    
363  c             The variables us, sh and rdn have to be computed from given  c             The variables us, sh and rdn have to be computed from
364  c             wind stresses inverting relationship for neutral drag coeff.  c             given wind stresses inverting relationship for neutral
365  c             cdn.  c             drag coeff. cdn.
366  c             The inversion is based on linear and quadratic form of  c             The inversion is based on linear and quadratic form of
367  c             cdn(umps); ustar can be directly computed from stress;  c             cdn(umps); ustar can be directly computed from stress;
368    
369                ustar = sqrt(ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +                ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +
370       &                     vstress(i,j,bi,bj)*vstress(i,j,bi,bj))/       &                vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
371       &                    atmrho                if ( ustmp .ne. 0. _d 0 ) then
372                cw = ustress(i,j,bi,bj)/ustar                  ustar = sqrt(ustmp/atmrho)
373                sw = ustress(i,j,bi,bj)/ustar                  cw = ustress(i,j,bi,bj)/sqrt(ustmp)
374                    sw = vstress(i,j,bi,bj)/sqrt(ustmp)
375                  else
376                     ustar = 0. _d 0
377                     cw    = 0. _d 0
378                     sw    = 0. _d 0
379                  endif
380    
381                if ( ustar .eq. 0. _d 0 ) then                if ( ustar .eq. 0. _d 0 ) then
382                  us = 0. _d 0                  us = 0. _d 0
# Line 362  c             cdn(umps); ustar can be di Line 404  c             cdn(umps); ustar can be di
404  #endif /* ALLOW_ATM_WIND */  #endif /* ALLOW_ATM_WIND */
405    
406  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
407    
408  c             Initial guess: z/l=0.0; hu=ht=hq=z  c             Initial guess: z/l=0.0; hu=ht=hq=z
409  c             Iterations:    converge on z/l and hence the fluxes.  c             Iterations:    converge on z/l and hence the fluxes.
410  c             t0     : virtual temperature (K)  c             t0     : virtual temperature (K)
# Line 371  c             deltap : potential tempera Line 414  c             deltap : potential tempera
414                if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then                if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then
415                  t0     = atemp(i,j,bi,bj)*                  t0     = atemp(i,j,bi,bj)*
416       &                   (exf_one + humid_fac*aqh(i,j,bi,bj))       &                   (exf_one + humid_fac*aqh(i,j,bi,bj))
417                    ssttmp = theta(i,j,k,bi,bj)
418                  ssq    = saltsat*                  ssq    = saltsat*
419       &                   exf_BulkqSat(theta(i,j,1,bi,bj) + cen2kel)/       &                   exf_BulkqSat(ssttmp + cen2kel)/
420       &                   atmrho       &                   atmrho
421                  deltap = atemp(i,j,bi,bj)   + gamma_blk*ht -                  deltap = atemp(i,j,bi,bj)   + gamma_blk*ht -
422       &                   theta(i,j,1,bi,bj) - cen2kel       &                   ssttmp - cen2kel
423                  delq   = aqh(i,j,bi,bj) - ssq                  delq   = aqh(i,j,bi,bj) - ssq
424                  stable = exf_half + sign(exf_half, deltap)                  stable = exf_half + sign(exf_half, deltap)
425    #ifdef ALLOW_AUTODIFF_TAMC
426    CADJ STORE sh     = comlev1_exf_1, key = ikey_1
427    #endif
428                  rdn    = sqrt(exf_BulkCdn(sh))                  rdn    = sqrt(exf_BulkCdn(sh))
429                  ustar  = rdn*sh                  ustar  = rdn*sh
430                  tstar  = exf_BulkRhn(stable)*deltap                  tstar  = exf_BulkRhn(stable)*deltap
# Line 385  c             deltap : potential tempera Line 432  c             deltap : potential tempera
432    
433                  do iter = 1,niter_bulk                  do iter = 1,niter_bulk
434    
435    #ifdef ALLOW_AUTODIFF_TAMC
436                       ikey_2 = iter
437         &                  + niter_bulk*(i-1)
438         &                  + sNx*niter_bulk*(j-1)
439         &                  + sNx*niter_bulk*sNy*act1
440         &                  + sNx*niter_bulk*sNy*max1*act2
441         &                  + sNx*niter_bulk*sNy*max1*max2*act3
442         &                  + sNx*niter_bulk*sNy*max1*max2*max3*act4
443    #endif
444    
445    #ifdef ALLOW_AUTODIFF_TAMC
446    CADJ STORE rdn    = comlev1_exf_2, key = ikey_2
447    CADJ STORE ustar  = comlev1_exf_2, key = ikey_2
448    CADJ STORE qstar  = comlev1_exf_2, key = ikey_2
449    CADJ STORE tstar  = comlev1_exf_2, key = ikey_2
450    CADJ STORE sh     = comlev1_exf_2, key = ikey_2
451    CADJ STORE us     = comlev1_exf_2, key = ikey_2
452    #endif
453    
454                    huol   = czol*(tstar/t0 +                    huol   = czol*(tstar/t0 +
455       &                     qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/       &                     qstar/(exf_one/humid_fac+aqh(i,j,bi,bj)))/
456       &                           ustar**2       &                           ustar**2
# Line 406  c                 Evaluate all stability Line 472  c                 Evaluate all stability
472       &                     exf_two*log((exf_one + xsq)/exf_two)       &                     exf_two*log((exf_one + xsq)/exf_two)
473    
474  c                 Shift wind speed using old coefficient  c                 Shift wind speed using old coefficient
475  c                 rd   = rdn/(exf_one + rdn/karman*(log(hu/zref) - psimh) )  ccc                  rd   = rdn/(exf_one + rdn/karman*
476    ccc     &                 (log(hu/zref) - psimh) )
477                    rd   = rdn/(exf_one - rdn/karman*psimh )                    rd   = rdn/(exf_one - rdn/karman*psimh )
478                    uzn  = max(sh*rd/rdn, umin)                    shn  = sh*rd/rdn
479                      uzn  = max(shn, umin)
480    
481  c                 Update the transfer coefficients at 10 meters  c                 Update the transfer coefficients at 10 meters
482  c                 and neutral stability.  c                 and neutral stability.
483    
484                    rdn = sqrt(exf_BulkCdn(uzn))                    rdn = sqrt(exf_BulkCdn(uzn))
485    
486  c                 Shift all coefficients to the measurement height  c                 Shift all coefficients to the measurement height
# Line 428  c                 coefficients. Line 497  c                 coefficients.
497                    ustar = rd*sh                      ustar = rd*sh  
498                    qstar = re*delq                    qstar = re*delq
499                    tstar = rh*deltap                    tstar = rh*deltap
500                      tau   = atmrho*ustar**2
501                      tau   = tau*us/sh
502    
503                  enddo                  enddo
504    
505                  tau                = atmrho*ustar**2    #ifdef ALLOW_AUTODIFF_TAMC
506                  tau                = tau*us/sh  CADJ STORE ustar  = comlev1_exf_1, key = ikey_1
507    CADJ STORE qstar  = comlev1_exf_1, key = ikey_1
508    CADJ STORE tstar  = comlev1_exf_1, key = ikey_1
509    CADJ STORE tau    = comlev1_exf_1, key = ikey_1
510    CADJ STORE cw     = comlev1_exf_1, key = ikey_1
511    CADJ STORE sw     = comlev1_exf_1, key = ikey_1
512    #endif
513    
514                  hs(i,j,bi,bj)      = atmcp*tau*tstar/ustar                  hs(i,j,bi,bj)      = atmcp*tau*tstar/ustar
515                  hl(i,j,bi,bj)      = flamb*tau*qstar/ustar                  hl(i,j,bi,bj)      = flamb*tau*qstar/ustar
516    
# Line 464  ce     &                             fla Line 542  ce     &                             fla
542        enddo        enddo
543    
544  c     Add all contributions.  c     Add all contributions.
       k = 1  
545        do bj = mybylo(mythid),mybyhi(mythid)        do bj = mybylo(mythid),mybyhi(mythid)
546          do bi = mybxlo(mythid),mybxhi(mythid)          do bi = mybxlo(mythid),mybxhi(mythid)
547            do j = 1,sny            do j = 1,sny
# Line 479  c             Net surface heat flux. Line 556  c             Net surface heat flux.
556                hfl = hfl + swflux(i,j,bi,bj)                hfl = hfl + swflux(i,j,bi,bj)
557  #endif /* ALLOW_KPP undef */  #endif /* ALLOW_KPP undef */
558  c             Heat flux:  c             Heat flux:
559                hflux(i,j,bi,bj) = hfl*maskc(i,j,k,bi,bj)                hflux(i,j,bi,bj) = hfl*maskc(i,j,1,bi,bj)
560  c             Salt flux from Precipitation and Evaporation.  c             Salt flux from Precipitation and Evaporation.
561                sflux(i,j,bi,bj) = precip(i,j,bi,bj) - evap(i,j,bi,bj)                sflux(i,j,bi,bj) = precip(i,j,bi,bj) - evap(i,j,bi,bj)
562  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
563    
564  #else  #else
565                hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j,k,bi,bj)                hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)
566                sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,k,bi,bj)                sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskc(i,j,1,bi,bj)
567  #endif /* ALLOW_BULKFORMULAE */  #endif /* ALLOW_BULKFORMULAE */
568              enddo              enddo
569            enddo            enddo
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22