/[MITgcm]/MITgcm/pkg/bulk_force/bulkf_formula_lanl.F
ViewVC logotype

Diff of /MITgcm/pkg/bulk_force/bulkf_formula_lanl.F

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

revision 1.7 by jmc, Sun Jan 22 15:51:35 2006 UTC revision 1.8 by jmc, Thu May 25 17:30:54 2006 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "BULK_FORCE_OPTIONS.h"  #include "BULK_FORCE_OPTIONS.h"
5    
6        subroutine bulkf_formula_lanl(  CBOP
7    C     !ROUTINE: BULKF_FORMULA_LANL
8    C     !INTERFACE:
9          SUBROUTINE BULKF_FORMULA_LANL(
10       I                           uw, vw, us, Ta, Qa, nc, tsf_in,       I                           uw, vw, us, Ta, Qa, nc, tsf_in,
11       O                           flwupa, flha, fsha, df0dT,       O                           flwupa, flha, fsha, df0dT,
12       O                           ust, vst, evp, ssq, dEvdT,       O                           ust, vst, evp, ssq, dEvdT,
13       I                           iceornot, myThid )       I                           iceornot, myThid )
14    
15  c swd -- bulkf formula used in bulkf and ice pkgs  C     !DESCRIPTION: \bv
16  c taken from exf package  C     *==========================================================*
17    C     | SUBROUTINE BULKF_FORMULA_LANL
18    c     | o Calculate bulk formula fluxes over open ocean or seaice
19    C     *==========================================================*
20    C     \ev
21    C swd -- bulkf formula used in bulkf and ice pkgs
22    C        taken from exf package
23    C
24    C     wind stress = (ust,vst) = rhoA * Cd * Ws * (del.u,del.v)
25    C     Sensib Heat flux = fsha = rhoA * Ch * Ws * del.T * CpAir
26    C     Latent Heat flux = flha = rhoA * Ce * Ws * del.Q * Lvap
27    C                      = -Evap * Lvap
28    C   with Ws = wind speed = sqrt(del.u^2 +del.v^2) ;
29    C        del.T = Tair - Tsurf ; del.Q = Qair - Qsurf
30    C        Cd,Ch,Ce = transfer coefficient for momentum, sensible
31    C                    & latent heat flux [no units]
32    C     *==========================================================*
33    
34    C     !USES:
35         IMPLICIT NONE         IMPLICIT NONE
36  c  C     === Global variables ===
37  #include "EEPARAMS.h"  #include "EEPARAMS.h"
38  #include "SIZE.h"  #include "SIZE.h"
39  #include "PARAMS.h"  #include "PARAMS.h"
40  #include "BULKF_PARAMS.h"  #include "BULKF_PARAMS.h"
41    
42  c  C     !INPUT/OUTPUT PARAMETERS:
43  c calculate bulk formula fluxes over open ocean or seaice  C     input:
44  c        _RL uw                 ! zonal wind speed (at grid center) [m/s]
45  c input        _RL vw                 ! meridional wind speed (at grid center) [m/s]
46        _RL us                 ! wind speed        _RL us                 ! wind speed        [m/s]   at height hu
47        _RL uw                 ! zonal wind speed (at grid center)        _RL Ta                 ! air temperature   [K]     at height ht
48        _RL vw                 ! meridional wind speed (at grid center)        _RL Qa                 ! specific humidity [kg/kg] at heigth ht
       _RL Ta                 ! air temperature at ht  
       _RL Qa                 ! specific humidity at ht  
       _RL tsf_in             ! surface temperature (either ice or open water)  
49        _RL nc                 ! fraction cloud cover        _RL nc                 ! fraction cloud cover
50        integer iceornot       ! 0=open water, 1=ice cover        _RL tsf_in             ! sea-ice or sea surface temperature [oC]
51        integer myThid         ! my Thread Id number        INTEGER iceornot       ! 0=open water, 1=sea-ice, 2=sea-ice with snow
52  c output        INTEGER myThid         ! my Thread Id number
53        _RL flwupa             ! upward long wave radiation (>0 upward)  C     output:
54        _RL flha               ! latent heat flux         (>0 downward)        _RL flwupa             ! upward long wave radiation (>0 upward) [W/m2]
55        _RL fsha               ! sensible heat flux       (>0 downward)        _RL flha               ! latent heat flux         (>0 downward) [W/m2]
56        _RL df0dT              ! derivative of heat flux with respect to Tsf        _RL fsha               ! sensible heat flux       (>0 downward) [W/m2]
57        _RL ust                ! zonal wind stress (at grid center)        _RL df0dT              ! derivative of heat flux with respect to Tsf [W/m2/K]
58        _RL vst                ! meridional wind stress (at grid center)        _RL ust                ! zonal wind stress (at grid center)     [N/m2]
59          _RL vst                ! meridional wind stress (at grid center)[N/m2]
60        _RL evp                ! evaporation rate (over open water) [kg/m2/s]        _RL evp                ! evaporation rate (over open water) [kg/m2/s]
61        _RL ssq                ! surface specific humidity (kg/kg)        _RL ssq                ! surface specific humidity          [kg/kg]
62        _RL dEvdT              ! derivative of evap. with respect to Tsf [kg/m2/s/K]        _RL dEvdT              ! derivative of evap. with respect to Tsf [kg/m2/s/K]
63    CEOP
64    
65  #ifdef ALLOW_BULK_FORCE  #ifdef ALLOW_BULK_FORCE
66    
67  c local variables    C     == Local variables ==
       _RL tsf                ! surface temperature in K  
       _RL ht                 ! reference height (m)  
       _RL hq                 ! reference height for humidity (m)  
       _RL hu                 ! reference height for wind speed (m)  
       _RL zref               ! reference height  
       _RL czol  
       _RL usm                ! wind speed limited  
       _RL t0                 ! virtual temperature (K)  
       _RL deltap             ! potential temperature diff (K)  
       _RL delq               ! specific humidity difference  
       _RL stable  
       _RL rdn   ,ren, rhn  
       _RL ustar  
       _RL tstar  
       _RL qstar  
       _RL huol  
       _RL xsq  
       _RL x  
       _RL re  
       _RL rh  
       _RL tau  
       _RL psimh  
       _RL psixh  
       _RL rd  
       _RL aln  
       _RL cdalton  
68        _RL dflhdT             ! derivative of latent heat with respect to T        _RL dflhdT             ! derivative of latent heat with respect to T
69        _RL dfshdT             ! derivative of sensible heat with respect to T        _RL dfshdT             ! derivative of sensible heat with respect to T
70        _RL dflwupdT           ! derivative of long wave with respect to T        _RL dflwupdT           ! derivative of long wave with respect to T
71    
72          _RL tsf                ! surface temperature [K]
73          _RL ht                 ! height for air temperature [m]
74    c     _RL hq                 ! height for humidity [m]
75          _RL hu                 ! height for wind speed [m]
76    c     _RL zref               ! reference height [m]
77          _RL usm                ! wind speed limited [m/s]
78    c     _RL umin               ! minimum wind speed used for drag-coeff [m/s]
79          _RL lath               ! latent heat of vaporization or sublimation
80          _RL t0                 ! virtual temperature [K]
81          _RL deltap             ! potential temperature diff [K]
82          _RL delq               ! specific humidity difference [kg/kg]
83          _RL ustar              ! friction velocity [m/s]
84          _RL tstar              ! temperature scale [K]
85          _RL qstar              ! humidity scale  [kg/kg]
86          _RL rd                 ! = sqrt(Cd)          [-]
87          _RL re                 ! = Ce / sqrt(Cd)     [-]
88          _RL rh                 ! = Ch / sqrt(Cd)     [-]
89          _RL rdn, ren, rhn      ! initial (neutral) values of rd, re, rh
90          _RL stable             ! = 1 if stable ; = 0 if unstable
91          _RL huol               ! stability parameter [-]
92          _RL x                  ! stability function  [-]
93          _RL xsq                ! = x^2               [-]
94          _RL psimh              ! momentum stability function
95          _RL psixh              ! latent & sensib. stability function
96          _RL czol               ! = zref*Karman_cst*gravity
97          _RL aln                ! = log(ht/zref)
98    c     _RL cdalton            ! coeff to evaluate Dalton Number
99  c     _RL mixratio  c     _RL mixratio
100  c     _RL ea  c     _RL ea
101        _RL psim_fac  c     _RL psim_fac
102        _RL umin        _RL tau                ! surface stress  coef = ?
103        _RL lath        _RL csha               ! sensib.heat flx coef = rhoA * Ws * Ch * CpAir
104        _RL csha        _RL clha               ! latent heat flx coef = rhoA * Ws * Ce * Lvap
       _RL clha  
105        _RL zice        _RL zice
106        _RL ssq0, ssq1, ssq2   ! constant used in surface specific humidity        _RL ssq0, ssq1, ssq2   ! constant used in surface specific humidity
107          _RL p0                 ! reference sea-level atmospheric pressure [mb]
108        _RL bulkf_Cdn          ! drag coefficient        _RL bulkf_Cdn          ! drag coefficient
109        integer niter_bulk, iter        INTEGER niter_bulk, iter
110    
111  c --- external functions ---  C     == external Functions
112  c     _RL       exf_BulkCdn  c     _RL       exf_BulkCdn
113  c     external  exf_BulkCdn  c     external  exf_BulkCdn
114  c     _RL       exf_BulkqSat  c     _RL       exf_BulkqSat
# Line 95  c     external  exf_BulkqSat Line 116  c     external  exf_BulkqSat
116  c     _RL       exf_BulkRhn  c     _RL       exf_BulkRhn
117  c     external  exf_BulkRhn  c     external  exf_BulkRhn
118    
119        DATA   ssq0,           ssq1,           ssq2        DATA   ssq0,           ssq1,           ssq2
120       &     / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /       &     / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /
121          DATA   p0 / 1013. _d 0 /
122    
123  cQQQQQQQQQQQQQQQQQQQQq  C--- Compute turbulent surface fluxes
 c -- compute turbulent surface fluxes  
124                ht =  2. _d 0                ht =  2. _d 0
125                hq =  2. _d 0  c             hq =  2. _d 0
126                hu = 10. _d 0                hu = 10. _d 0
127                zref = 10. _d 0  c             zref = 10. _d 0
128                zice = 0.0005 _d 0                zice = 0.0005 _d 0
129                aln = log(ht/zref)                aln = log(ht/zref)
130                niter_bulk = 5                niter_bulk = 5
131                cdalton = 0.0346000 _d 0  c             cdalton = 0.0346000 _d 0
132                czol = zref*xkar*gravity                czol = zref*xkar*gravity
133                psim_fac=5. _d 0  c             psim_fac=5. _d 0
134                umin=1. _d 0  c             umin=1. _d 0
135  c  
136                lath=Lvap                lath=Lvap
137                if (iceornot.gt.0) lath=Lvap+Lfresh                if (iceornot.gt.0) lath=Lvap+Lfresh
138                Tsf=Tsf_in+Tf0kel                Tsf=Tsf_in+Tf0kel
139  c     Wind speed  C-   Wind speed
140                if (us.eq.0. _d 0) then                if (us.eq.0. _d 0) then
141                  us = sqrt(uw*uw + vw*vw)                  us = sqrt(uw*uw + vw*vw)
142                endif                endif
143                usm = max(us,umin)                usm = max(us,umin)
 cQQQ try to limit drag  
 cQQ           usm = min(usm,5. _d 0)  
144  c  c
145                t0     = Ta*(1. _d 0 + humid_fac*Qa)                t0     = Ta*(1. _d 0 + humid_fac*Qa)
 cQQ           ssq    = 0.622*6.11*exp(22.47*(1.d0-Tf0kel/tsf))/p0  
146  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
147    cQQ           ssq    = 0.622*6.11*exp(22.47*(1.d0-Tf0kel/tsf))/p0
148  c             ssq = 3.797915 _d 0*exp(  c             ssq = 3.797915 _d 0*exp(
149  c    &                lath*(7.93252 _d -6 - 2.166847 _d -3/Tsf)  c    &                lath*(7.93252 _d -6 - 2.166847 _d -3/Tsf)
150  c    &                               )/p0  c    &                               )/p0
151                ssq = ssq0*exp( lath*(ssq1-ssq2/Tsf) ) / p0                ssq = ssq0*exp( lath*(ssq1-ssq2/Tsf) ) / p0
152  cBB debugging  
 cBB             print*,'ice, ssq',  ssq  
 c  
153                deltap = ta  - tsf + gamma_blk*ht                deltap = ta  - tsf + gamma_blk*ht
154                delq   = Qa - ssq                delq   = Qa - ssq
155  c  
156  c initialize estimate exchange coefficients  C--  initialize estimate exchange coefficients
157                rdn=xkar/(log(zref/zice))                rdn=xkar/(log(zref/zice))
158                rhn=rdn                rhn=rdn
159                ren=rdn                ren=rdn
160  c calculate turbulent scales  C--  calculate turbulent scales
161                ustar=rdn*usm                ustar=rdn*usm
162                tstar=rhn*deltap                tstar=rhn*deltap
163                qstar=ren*delq                qstar=ren*delq
164  c  
165  c interation with psi-functions to find transfer coefficients  C--  iteration with psi-functions to find transfer coefficients
166                do iter=1,niter_bulk                do iter=1,niter_bulk
167                   huol   = czol/ustar**2 *(tstar/t0 +                   huol   = czol/ustar**2 *(tstar/t0 +
168       &                    qstar/(1. _d 0/humid_fac+Qa))       &                    qstar/(1. _d 0/humid_fac+Qa))
# Line 159  c interation with psi-functions to find Line 176  c interation with psi-functions to find
176       &                     2. _d 0*atan(x) + pi*.5 _d 0)       &                     2. _d 0*atan(x) + pi*.5 _d 0)
177                   psixh  = -5. _d 0*huol*stable + (1. _d 0-stable)*                   psixh  = -5. _d 0*huol*stable + (1. _d 0-stable)*
178       &                     (2. _d 0*log(5. _d -1*(1. _d 0+xsq)))       &                     (2. _d 0*log(5. _d -1*(1. _d 0+xsq)))
179      
180  c Update the transfer coefficients  C--  Update the transfer coefficients
     
181                   rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar)                   rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar)
182                   rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar)                   rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar)
183                   re = rh                   re = rh
184  c  Update ustar, tstar, qstar using updated, shifted coefficients.  C--  Update ustar, tstar, qstar using updated, shifted coefficients.
185                   ustar = rd*usm                   ustar = rd*usm
186                   qstar = re*delq                   qstar = re*delq
187                   tstar = rh*deltap                   tstar = rh*deltap
188                enddo                enddo
189  c  
190                  tau   = rhoa*ustar**2                tau   = rhoa*ustar**2
191                  tau   = tau*us/usm                tau   = tau*us/usm
192                  csha   = rhoa*cpair*us*rh*rd                  csha  = rhoa*cpair*us*rh*rd
193                  clha   = rhoa*lath*us*re*rd                  clha  = rhoa*lath*us*re*rd
194  c  
195                  fsha  = csha*deltap                fsha  = csha*deltap
196                  flha  = clha*delq                flha  = clha*delq
197                  evp   = -flha/lath                evp   = -flha/lath
198  c  
199  c up long wave radiation  C--  Upward long wave radiation
200  cQQ           mixratio=Qa/(1-Qa)  cQQ           mixratio=Qa/(1-Qa)
201  cQQ           ea=p0*mixratio/(0.62197+mixratio)  cQQ           ea=p0*mixratio/(0.62197+mixratio)
202  cQQ           flwupa=-0.985*stefan*tsf**4  cQQ           flwupa=-0.985*stefan*tsf**4
# Line 189  cQQ  &                  *(1-0.6*nc**2) Line 205  cQQ  &                  *(1-0.6*nc**2)
205                if (iceornot.eq.0) then                if (iceornot.eq.0) then
206                 flwupa=ocean_emissivity*stefan*tsf**4                 flwupa=ocean_emissivity*stefan*tsf**4
207                 dflwupdT=4. _d 0*ocean_emissivity*stefan*tsf**3                 dflwupdT=4. _d 0*ocean_emissivity*stefan*tsf**3
208                else                elseif (iceornot.eq.2) then
                if (iceornot.eq.2) then  
209                  flwupa=snow_emissivity*stefan*tsf**4                  flwupa=snow_emissivity*stefan*tsf**4
210                  dflwupdT=4. _d 0*snow_emissivity*stefan*tsf**3                  dflwupdT=4. _d 0*snow_emissivity*stefan*tsf**3
211                 else                else
212                  flwupa=ice_emissivity*stefan*tsf**4                  flwupa=ice_emissivity*stefan*tsf**4
213                  dflwupdT=4. _d 0*ice_emissivity*stefan*tsf**3                  dflwupdT=4. _d 0*ice_emissivity*stefan*tsf**3
                endif  
214                endif                endif
215  cQQ           dflhdT = -clha*Tf0kel*ssq*22.47/(tsf**2)  cQQ           dflhdT = -clha*Tf0kel*ssq*22.47/(tsf**2)
216  c             dflhdT = -clha*Lath*ssq/(Rvap*tsf**2)  c             dflhdT = -clha*Lath*ssq/(Rvap*tsf**2)
# Line 209  cQQ  &                  *(0.39-0.05*sqrt Line 223  cQQ  &                  *(0.39-0.05*sqrt
223  cQQ  &                  *(1-0.6*nc**2)  cQQ  &                  *(1-0.6*nc**2)
224  c total derivative with respect to surface temperature  c total derivative with respect to surface temperature
225                df0dT=-dflwupdT+dfshdT+dflhdT                df0dT=-dflwupdT+dfshdT+dflhdT
226  c  
227  c wind stress at center points  C--  Wind stress at center points
228  c             if (.NOT.windread) then  C-   in-lining of function: exf_BulkCdn(umps) = cdrag_1/umps + cdrag_2 + cdrag_3*umps
 c                ust = rhoa*exf_BulkCdn(usm)*us*uw  
 c                vst = rhoa*exf_BulkCdn(usm)*us*vw  
 c             endif  
 C-  In-lining of function: exf_BulkCdn(umps) = cdrag_1/umps + cdrag_2 + cdrag_3*umps  
229                bulkf_Cdn = cdrag_1/usm + cdrag_2 + cdrag_3*usm                bulkf_Cdn = cdrag_1/usm + cdrag_2 + cdrag_3*usm
230                ust = rhoa*bulkf_Cdn*us*uw                ust = rhoa*bulkf_Cdn*us*uw
231                vst = rhoa*bulkf_Cdn*us*vw                vst = rhoa*bulkf_Cdn*us*vw
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22