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

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

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

revision 1.11 by heimbach, Tue Dec 13 19:46:46 2005 UTC revision 1.12 by heimbach, Thu May 25 18:32:55 2006 UTC
# Line 2  c $Header$ Line 2  c $Header$
2    
3  #include "EXF_OPTIONS.h"  #include "EXF_OPTIONS.h"
4    
5        subroutine exf_bulkformulae(mycurrenttime, mycurrentiter, mythid)        subroutine exf_bulkformulae(mytime, myiter, mythid)
6    
7  c     ==================================================================  c     ==================================================================
8  c     SUBROUTINE exf_bulkformulae  c     SUBROUTINE exf_bulkformulae
# Line 128  c     == global variables == Line 128  c     == global variables ==
128  c     == routine arguments ==  c     == routine arguments ==
129    
130        integer mythid        integer mythid
131        integer mycurrentiter        integer myiter
132        _RL     mycurrenttime        _RL     mytime
133    
134  #ifdef ALLOW_BULKFORMULAE  #ifdef ALLOW_BULKFORMULAE
135    
# Line 171  c     == local variables == Line 171  c     == local variables ==
171  #endif  #endif
172  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
173    
       _RL     ustmp  
       _RL     us  
       _RL     cw  
       _RL     sw  
       _RL     sh  
174        _RL     hfl        _RL     hfl
175        _RL     tmpbulk        _RL     tmpbulk
176    
# Line 191  c     == external functions == Line 186  c     == external functions ==
186        _RL       exf_BulkRhn        _RL       exf_BulkRhn
187        external  exf_BulkRhn        external  exf_BulkRhn
188    
 #ifndef ALLOW_ATM_WIND  
       _RL       TMP1  
       _RL       TMP2  
       _RL       TMP3  
       _RL       TMP4  
       _RL       TMP5  
 #endif  
   
189  c     == end of interface ==  c     == end of interface ==
190    
191  cph   This statement cannot be a PARAMETER statement in the header,  cph   This statement cannot be a PARAMETER statement in the header,
192  cph   but must come here; it's not fortran77 standard  cph   but must come here; it is not fortran77 standard
193        aln = log(ht/zref)        aln = log(ht/zref)
194    
195  c--   Use atmospheric state to compute surface fluxes.  c--   Use atmospheric state to compute surface fluxes.
# Line 217  CHPF$ INDEPENDENT Line 204  CHPF$ INDEPENDENT
204  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
205  CHPF$  INDEPENDENT  CHPF$  INDEPENDENT
206  #endif  #endif
207          do bi = mybxlo(mythid),mybxhi(mythid)         do bi = mybxlo(mythid),mybxhi(mythid)
208            k = 1
209            k = 1          do j = 1,sny
210             do i = 1,snx
           do j = 1,sny  
             do i = 1,snx  
211    
212  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
213                 act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
214                 max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
215                 act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
216                 max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
217                 act3 = myThid - 1            act3 = myThid - 1
218                 max3 = nTx*nTy            max3 = nTx*nTy
219                 act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
220    
221                 ikey_1 = i            ikey_1 = i
222       &              + sNx*(j-1)       &         + sNx*(j-1)
223       &              + sNx*sNy*act1       &         + sNx*sNy*act1
224       &              + sNx*sNy*max1*act2       &         + sNx*sNy*max1*act2
225       &              + sNx*sNy*max1*max2*act3       &         + sNx*sNy*max1*max2*act3
226       &              + sNx*sNy*max1*max2*max3*act4       &         + sNx*sNy*max1*max2*max3*act4
 #endif  
   
 #ifdef ALLOW_DOWNWARD_RADIATION  
 c--   Compute net from downward and downward from net longwave and  
 c     shortwave radiation, if needed.  
 c     lwflux = Stefan-Boltzmann constant * emissivity * SST - lwdown  
 c     swflux = - ( 1 - albedo ) * swdown  
   
 #ifdef ALLOW_ATM_TEMP  
       if ( lwfluxfile .EQ. ' ' .AND. lwdownfile .NE. ' ' )  
      &     lwflux(i,j,bi,bj) = 5.5 _d -08 *  
      &              ((theta(i,j,k,bi,bj)+cen2kel)**4)  
      &              - lwdown(i,j,bi,bj)  
       if ( lwfluxfile .NE. ' ' .AND. lwdownfile .EQ. ' ' )  
      &     lwdown(i,j,bi,bj) = 5.5 _d -08 *  
      &              ((theta(i,j,k,bi,bj)+cen2kel)**4)  
      &              - lwflux(i,j,bi,bj)  
 #endif  
   
 #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)  
       if ( swfluxfile .EQ. ' ' .AND. swdownfile .NE. ' ' )  
      &     swflux(i,j,bi,bj) = -(1.0-exf_albedo) * swdown(i,j,bi,bj)  
       if ( swfluxfile .NE. ' ' .AND. swdownfile .EQ. ' ' )  
      &     swdown(i,j,bi,bj) = -swflux(i,j,bi,bj) / (1.0-exf_albedo)  
227  #endif  #endif
228    
 #endif /* ALLOW_DOWNWARD_RADIATION */  
   
229  c--   Compute the turbulent surface fluxes.  c--   Compute the turbulent surface fluxes.
230    
 #ifdef ALLOW_ATM_WIND  
 c             Wind speed and direction.  
               ustmp = uwind(i,j,bi,bj)*uwind(i,j,bi,bj) +  
      &                vwind(i,j,bi,bj)*vwind(i,j,bi,bj)  
               if ( ustmp .ne. 0. _d 0 ) then  
                 us = sqrt(ustmp)  
                 cw = uwind(i,j,bi,bj)/us  
                 sw = vwind(i,j,bi,bj)/us  
               else  
                 us = 0. _d 0  
                 cw = 0. _d 0  
                 sw = 0. _d 0  
               endif  
               sh = max(us,umin)  
 #else  /* ifndef ALLOW_ATM_WIND */  
 #ifdef ALLOW_ATM_TEMP  
   
 c             The variables us, sh and rdn have to be computed from  
 c             given wind stresses inverting relationship for neutral  
 c             drag coeff. cdn.  
 c             The inversion is based on linear and quadratic form of  
 c             cdn(umps); ustar can be directly computed from stress;  
   
               ustmp = ustress(i,j,bi,bj)*ustress(i,j,bi,bj) +  
      &                vstress(i,j,bi,bj)*vstress(i,j,bi,bj)  
               if ( ustmp .ne. 0. _d 0 ) then  
                 ustar = sqrt(ustmp/atmrho)  
                 cw = ustress(i,j,bi,bj)/sqrt(ustmp)  
                 sw = vstress(i,j,bi,bj)/sqrt(ustmp)  
               else  
                  ustar = 0. _d 0  
                  cw    = 0. _d 0  
                  sw    = 0. _d 0  
               endif  
   
               if ( ustar .eq. 0. _d 0 ) then  
                 us = 0. _d 0  
               else if ( ustar .lt. ustofu11 ) then  
                 tmp1 = -cquadrag_2/cquadrag_1/2  
                 tmp2 = sqrt(tmp1*tmp1 + ustar*ustar/cquadrag_1)  
                 us   = sqrt(tmp1 + tmp2)  
               else  
                 tmp3 = clindrag_2/clindrag_1/3  
                 tmp4 = ustar*ustar/clindrag_1/2 - tmp3**3  
                 tmp5 = sqrt(ustar*ustar/clindrag_1*  
      &                      (ustar*ustar/clindrag_1/4 - tmp3**3))  
                 us   = (tmp4 + tmp5)**(1/3) +  
      &                 tmp3**2 * (tmp4 + tmp5)**(-1/3) - tmp3  
               endif  
   
               if ( us .ne. 0 ) then  
                  rdn = ustar/us  
               else  
                  rdn = 0. _d 0  
               end if  
   
               sh    = max(us,umin)  
 #endif /* ALLOW_ATM_TEMP */  
 #endif /* ifndef ALLOW_ATM_WIND */  
   
231  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
232    
233  c             Initial guess: z/l=0.0; hu=ht=hq=z  c             Initial guess: z/l=0.0; hu=ht=hq=z
# Line 337  c             t0     : virtual temperatu Line 236  c             t0     : virtual temperatu
236  c             ssq    : sea surface humidity (kg/kg)  c             ssq    : sea surface humidity (kg/kg)
237  c             deltap : potential temperature diff (K)  c             deltap : potential temperature diff (K)
238    
239                if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then            if ( atemp(i,j,bi,bj) .ne. 0. _d 0 ) then
240                  t0     = atemp(i,j,bi,bj)*               t0     = atemp(i,j,bi,bj)*
241       &                   (exf_one + humid_fac*aqh(i,j,bi,bj))       &            (exf_one + humid_fac*aqh(i,j,bi,bj))
242                  ssttmp = theta(i,j,k,bi,bj)               ssttmp = theta(i,j,k,bi,bj)
243                  tmpbulk= exf_BulkqSat(ssttmp + cen2kel)               tmpbulk= exf_BulkqSat(ssttmp + cen2kel)
244                  ssq    = saltsat*tmpbulk/atmrho               ssq    = saltsat*tmpbulk/atmrho
245                  deltap = atemp(i,j,bi,bj)   + gamma_blk*ht -               deltap = atemp(i,j,bi,bj)   + gamma_blk*ht -
246       &                   ssttmp - cen2kel       &            ssttmp - cen2kel
247                  delq   = aqh(i,j,bi,bj) - ssq               delq   = aqh(i,j,bi,bj) - ssq
248                  stable = exf_half + sign(exf_half, deltap)               stable = exf_half + sign(exf_half, deltap)
249  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
250  CADJ STORE sh     = comlev1_exf_1, key = ikey_1  CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1
251  #endif  #endif
252                  tmpbulk= exf_BulkCdn(sh)               tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))
253                  rdn    = sqrt(tmpbulk)               rdn    = sqrt(tmpbulk)
254                  ustar  = rdn*sh               ustar  = rdn*sh(i,j,bi,bj)
255                  tmpbulk= exf_BulkRhn(stable)               tmpbulk= exf_BulkRhn(stable)
256                  tstar  = tmpbulk*deltap               tstar  = tmpbulk*deltap
257                  qstar  = cdalton*delq               qstar  = cdalton*delq
258    
259                  do iter = 1,niter_bulk               do iter = 1,niter_bulk
260    
261  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
262                     ikey_2 = iter                     ikey_2 = iter
# Line 372  CADJ STORE rdn    = comlev1_exf_2, key = Line 271  CADJ STORE rdn    = comlev1_exf_2, key =
271  CADJ STORE ustar  = comlev1_exf_2, key = ikey_2  CADJ STORE ustar  = comlev1_exf_2, key = ikey_2
272  CADJ STORE qstar  = comlev1_exf_2, key = ikey_2  CADJ STORE qstar  = comlev1_exf_2, key = ikey_2
273  CADJ STORE tstar  = comlev1_exf_2, key = ikey_2  CADJ STORE tstar  = comlev1_exf_2, key = ikey_2
274  CADJ STORE sh     = comlev1_exf_2, key = ikey_2  CADJ STORE sh(i,j,bi,bj)     = comlev1_exf_2, key = ikey_2
275  CADJ STORE us     = comlev1_exf_2, key = ikey_2  CADJ STORE us(i,j,bi,bj)     = comlev1_exf_2, key = ikey_2
276  #endif  #endif
277    
278                    huol   = czol*(tstar/t0 +                    huol   = czol*(tstar/t0 +
# Line 400  c                 Shift wind speed using Line 299  c                 Shift wind speed using
299  ccc                  rd   = rdn/(exf_one + rdn/karman*  ccc                  rd   = rdn/(exf_one + rdn/karman*
300  ccc     &                 (log(hu/zref) - psimh) )  ccc     &                 (log(hu/zref) - psimh) )
301                    rd     = rdn/(exf_one - rdn/karman*psimh )                    rd     = rdn/(exf_one - rdn/karman*psimh )
302                    shn    = sh*rd/rdn                    shn    = sh(i,j,bi,bj)*rd/rdn
303                    uzn    = max(shn, umin)                    uzn    = max(shn, umin)
304    
305  c                 Update the transfer coefficients at 10 meters  c                 Update the transfer coefficients at 10 meters
# Line 421  c                 rd = rdn/(exf_one + rd Line 320  c                 rd = rdn/(exf_one + rd
320    
321  c                 Update ustar, tstar, qstar using updated, shifted  c                 Update ustar, tstar, qstar using updated, shifted
322  c                 coefficients.  c                 coefficients.
323                    ustar  = rd*sh                      ustar  = rd*sh(i,j,bi,bj)
324                    qstar  = re*delq                    qstar  = re*delq
325                    tstar  = rh*deltap                    tstar  = rh*deltap
326                    tau    = atmrho*ustar**2                    tau    = atmrho*ustar**2
327                    tau    = tau*us/sh                    tau    = tau*us(i,j,bi,bj)/sh(i,j,bi,bj)
328    
329                  enddo                  enddo
330    
# Line 434  CADJ STORE ustar  = comlev1_exf_1, key = Line 333  CADJ STORE ustar  = comlev1_exf_1, key =
333  CADJ STORE qstar  = comlev1_exf_1, key = ikey_1  CADJ STORE qstar  = comlev1_exf_1, key = ikey_1
334  CADJ STORE tstar  = comlev1_exf_1, key = ikey_1  CADJ STORE tstar  = comlev1_exf_1, key = ikey_1
335  CADJ STORE tau    = comlev1_exf_1, key = ikey_1  CADJ STORE tau    = comlev1_exf_1, key = ikey_1
336  CADJ STORE cw     = comlev1_exf_1, key = ikey_1  CADJ STORE cw(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1
337  CADJ STORE sw     = comlev1_exf_1, key = ikey_1  CADJ STORE sw(i,j,bi,bj)     = comlev1_exf_1, key = ikey_1
338  #endif  #endif
339    
340                  hs(i,j,bi,bj)      = atmcp*tau*tstar/ustar                  hs(i,j,bi,bj)      = atmcp*tau*tstar/ustar
# Line 445  cdm             evap(i,j,bi,bj)    = tau Line 344  cdm             evap(i,j,bi,bj)    = tau
344  cdm !!! need to change sign and to convert from kg/m^2/s to m/s !!!  cdm !!! need to change sign and to convert from kg/m^2/s to m/s !!!
345                  evap(i,j,bi,bj)    = -recip_rhonil*tau*qstar/ustar                  evap(i,j,bi,bj)    = -recip_rhonil*tau*qstar/ustar
346  #endif  #endif
347                  ustress(i,j,bi,bj) = tau*cw                  ustress(i,j,bi,bj) = tau*cw(i,j,bi,bj)
348                  vstress(i,j,bi,bj) = tau*sw                  vstress(i,j,bi,bj) = tau*sw(i,j,bi,bj)
349                else                else
350                  ustress(i,j,bi,bj) = 0. _d 0                  ustress(i,j,bi,bj) = 0. _d 0
351                  vstress(i,j,bi,bj) = 0. _d 0                  vstress(i,j,bi,bj) = 0. _d 0
# Line 457  cdm !!! need to change sign and to conve Line 356  cdm !!! need to change sign and to conve
356    
357  #else  /* ifndef ALLOW_ATM_TEMP */  #else  /* ifndef ALLOW_ATM_TEMP */
358  #ifdef ALLOW_ATM_WIND  #ifdef ALLOW_ATM_WIND
359                tmpbulk= exf_BulkCdn(sh)                tmpbulk= exf_BulkCdn(sh(i,j,bi,bj))
360                ustress(i,j,bi,bj) = atmrho*tmpbulk*us*                ustress(i,j,bi,bj) = atmrho*tmpbulk*us(i,j,bi,bj)*
361       &                             uwind(i,j,bi,bj)       &                             uwind(i,j,bi,bj)
362                vstress(i,j,bi,bj) = atmrho*tmpbulk*us*                vstress(i,j,bi,bj) = atmrho*tmpbulk*us(i,j,bi,bj)*
363       &                             vwind(i,j,bi,bj)       &                             vwind(i,j,bi,bj)
364  #endif  #endif
365  #endif /* ifndef ALLOW_ATM_TEMP */  #endif /* ifndef ALLOW_ATM_TEMP */
366              enddo              enddo
367            enddo            enddo
368          enddo          enddo
       enddo  
   
 c     Add all contributions.  
       do bj = mybylo(mythid),mybyhi(mythid)  
         do bi = mybxlo(mythid),mybxhi(mythid)  
           do j = 1,sny  
             do i = 1,snx  
 c             Net surface heat flux.  
 #ifdef ALLOW_ATM_TEMP  
               hfl = 0. _d 0  
               hfl = hfl - hs(i,j,bi,bj)  
               hfl = hfl - hl(i,j,bi,bj)  
               hfl = hfl + lwflux(i,j,bi,bj)  
 #ifndef SHORTWAVE_HEATING  
               hfl = hfl + swflux(i,j,bi,bj)  
 #endif  
 c             Heat flux:  
               hflux(i,j,bi,bj) = hfl  
 c             Salt flux from Precipitation and Evaporation.  
               sflux(i,j,bi,bj) = evap(i,j,bi,bj) - precip(i,j,bi,bj)  
 #endif /* ALLOW_ATM_TEMP */  
   
             enddo  
           enddo  
         enddo  
369        enddo        enddo
370    
371  #endif /* ALLOW_BULKFORMULAE */  #endif /* ALLOW_BULKFORMULAE */
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22