/[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.3 by edhill, Thu Oct 9 04:19:19 2003 UTC revision 1.4 by jmc, Sun Nov 23 01:36:55 2003 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
 c swd -- bulkf formula used in bulkf and ice pkgs  
3    
 c taken from exf package  
4  #include "BULK_FORCE_OPTIONS.h"  #include "BULK_FORCE_OPTIONS.h"
5    
6        subroutine bulkf_formula_lanl(        subroutine bulkf_formula_lanl(
7       I                           uw, vw, us, Ta, Qa, nc, tsf_in,       I                           uw, vw, us, Ta, Qa, nc, tsf_in,
8       I                           flwupa, flha, fsha, df0dT,       O                           flwupa, flha, fsha, df0dT,
9       I                           ust, vst, evp, ssq, iceornot, windread       O                           ust, vst, evp, ssq,
10       &                           )       I                           iceornot, windread )
11    
12    c swd -- bulkf formula used in bulkf and ice pkgs
13    c taken from exf package
14    
15         IMPLICIT NONE         IMPLICIT NONE
16  c  c
17  #include "EEPARAMS.h"  #include "EEPARAMS.h"
18  #include "SIZE.h"  #include "SIZE.h"
19  #include "PARAMS.h"  #include "PARAMS.h"
20  #include "DYNVARS.h"  #include "BULKF_PARAMS.h"
 #include "GRID.h"  
 #include "FFIELDS.h"  
 #ifdef ALLOW_BULK_FORCE  
 #include "BULKF_ICE_CONSTANTS.h"  
 #include "BULKF.h"  
 #endif  
21    
22  c  c
23  c calculate bulk formula fluxes over open ocean or seaice  c calculate bulk formula fluxes over open ocean or seaice
# Line 37  c input Line 33  c input
33        integer iceornot       ! 0=open water, 1=ice cover        integer iceornot       ! 0=open water, 1=ice cover
34        logical windread       !        logical windread       !
35  c output  c output
36        _RL flwupa              ! upward long wave radiation        _RL flwupa             ! upward long wave radiation (>0 upward)
37        _RL flha                ! latent heat flux        _RL flha               ! latent heat flux         (>0 downward)
38        _RL fsha                ! sensible heat flux        _RL fsha               ! sensible heat flux       (>0 downward)
39        _RL df0dT              ! derivative of heat flux with respect to Tsf        _RL df0dT              ! derivative of heat flux with respect to Tsf
40        _RL ust                ! zonal wind stress (at grid center)        _RL ust                ! zonal wind stress (at grid center)
41        _RL vst                ! meridional wind stress (at grid center)        _RL vst                ! meridional wind stress (at grid center)
42        _RL evp               ! evaporation rate (over open water)        _RL evp                ! evaporation rate (over open water)
43        _RL ssq                ! surface specific humidity (kg/kg)        _RL ssq                ! surface specific humidity (kg/kg)
44  c  
45    #ifdef ALLOW_BULK_FORCE
46    
47  c local variables    c local variables  
48        _RL tsf                ! surface temperature in K        _RL tsf                ! surface temperature in K
49        _RL ht                 ! reference height (m)        _RL ht                 ! reference height (m)
# Line 90  c local variables Line 88  c local variables
88        _RL zice        _RL zice
89        integer niter_bulk, iter        integer niter_bulk, iter
90    
 #ifdef ALLOW_BULK_FORCE  
   
91  c --- external functions ---  c --- external functions ---
92        _RL       exf_BulkCdn        _RL       exf_BulkCdn
93        external  exf_BulkCdn        external  exf_BulkCdn
94        _RL       exf_BulkqSat  c     _RL       exf_BulkqSat
95        external  exf_BulkqSat  c     external  exf_BulkqSat
96        _RL       exf_BulkRhn  c     _RL       exf_BulkRhn
97        external  exf_BulkRhn  c     external  exf_BulkRhn
98    
99  cQQQQQQQQQQQQQQQQQQQQq  cQQQQQQQQQQQQQQQQQQQQq
100  c -- compute turbulent surface fluxes  c -- compute turbulent surface fluxes
101                ht =  2.d0                ht =  2. _d 0
102                hq =  2.d0                hq =  2. _d 0
103                hu = 10.d0                hu = 10. _d 0
104                zref = 10.d0                zref = 10. _d 0
105                zice = 0.0005                zice = 0.0005 _d 0
106                aln = log(ht/zref)                aln = log(ht/zref)
107                niter_bulk = 5                niter_bulk = 5
108                cdalton = 0.0346000                cdalton = 0.0346000 _d 0
109                czol = zref*xkar*gravity                czol = zref*xkar*gravity
110                psim_fac=5.d0                psim_fac=5. _d 0
111                umin=1.d0                umin=1. _d 0
112  c  c
113                lath=Lvap                lath=Lvap
114                if (iceornot.gt.0) lath=Lvap+Lfresh                if (iceornot.gt.0) lath=Lvap+Lfresh
115                Tsf=Tsf_in+Tf0kel                Tsf=Tsf_in+Tf0kel
116  c     Wind speed  c     Wind speed
117                if (us.eq.0) then                if (us.eq.0. _d 0) then
118                  us = sqrt(uw*uw + vw*vw)                  us = sqrt(uw*uw + vw*vw)
119                endif                endif
120                usm = max(us,umin)                usm = max(us,umin)
121  cQQQ try to limit drag  cQQQ try to limit drag
122  cQQ           usm = min(usm,5.d0)  cQQ           usm = min(usm,5. _d 0)
123  c  c
124                t0     = Ta*(1.d0 + humid_fac*Qa)                t0     = Ta*(1. _d 0 + humid_fac*Qa)
125  cQQ           ssq    = 0.622*6.11*exp(22.47*(1.d0-Tf0kel/tsf))/p0  cQQ           ssq    = 0.622*6.11*exp(22.47*(1.d0-Tf0kel/tsf))/p0
126           ssq    =3.797915*exp(lath*(7.93252d-6-2.166847d-3/Tsf))/p0  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
127                  ssq = 3.797915 _d 0*exp(
128         &                lath*(7.93252 _d -6 - 2.166847 _d -3/Tsf)
129         &                               )/p0
130  cBB debugging  cBB debugging
131  cBB             print*,'ice, ssq',  ssq  cBB             print*,'ice, ssq',  ssq
132  c  c
# Line 146  c Line 145  c
145  c interation with psi-functions to find transfer coefficients  c interation with psi-functions to find transfer coefficients
146                do iter=1,niter_bulk                do iter=1,niter_bulk
147                   huol   = czol/ustar**2 *(tstar/t0 +                   huol   = czol/ustar**2 *(tstar/t0 +
148       &                    qstar/(1.d0/humid_fac+Qa))       &                    qstar/(1. _d 0/humid_fac+Qa))
149                   huol   = sign( min(abs(huol),10.d0), huol)                   huol   = sign( min(abs(huol),10. _d 0), huol)
150                   stable = 5. _d -1 + sign(5. _d -1 , huol)                   stable = 5. _d -1 + sign(5. _d -1 , huol)
151                   xsq    = max(sqrt(abs(1.d0 - 16.d0*huol)),1.d0)                   xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*huol)),1. _d 0)
152                   x      = sqrt(xsq)                   x      = sqrt(xsq)
153                   psimh = -5.d0*huol*stable + (1.d0-stable)*                   psimh = -5. _d 0*huol*stable + (1. _d 0-stable)*
154       &                    (2.d0*log(5.e-1*(1.d0+x)) +       &                    (2. _d 0*log(5. _d -1*(1. _d 0+x)) +
155       &                    2.d0*log(5.e-1*(1.d0+xsq)) -       &                     2. _d 0*log(5. _d -1*(1. _d 0+xsq)) -
156       &                    2.d0*atan(x) + pi*.5d0)       &                     2. _d 0*atan(x) + pi*.5 _d 0)
157                   psixh  = -5.d0*huol*stable + (1.d0-stable)*                   psixh  = -5. _d 0*huol*stable + (1. _d 0-stable)*
158       &                     (2.d0*log(5.e-1*(1.d0+xsq)))       &                     (2. _d 0*log(5. _d -1*(1. _d 0+xsq)))
159        
160  c Update the transfer coefficients  c Update the transfer coefficients
161        
162                   rd = rdn/(1.d0 + rdn*(aln-psimh)/xkar)                   rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar)
163                   rh = rhn/(1.d0 + rhn*(aln-psixh)/xkar)                   rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar)
164                   re = rh                   re = rh
165  c  Update ustar, tstar, qstar using updated, shifted coefficients.  c  Update ustar, tstar, qstar using updated, shifted coefficients.
166                   ustar = rd*usm                   ustar = rd*usm
167                   qstar = re*delq                   qstar = re*delq
168                   tstar = rh*deltap                   tstar = rh*deltap
169                  enddo                enddo
170  c  c
171                  tau   = rhoa*ustar**2                  tau   = rhoa*ustar**2
172                  tau   = tau*us/usm                  tau   = tau*us/usm
# Line 176  c Line 175  c
175  c  c
176                  fsha  = csha*deltap                  fsha  = csha*deltap
177                  flha  = clha*delq                  flha  = clha*delq
178                  evp   = flha/lath/rhosw                  evp   = -flha/lath/rhofw
179  c  c
180  c up long wave radiation  c up long wave radiation
181  cQQ           mixratio=Qa/(1-Qa)  cQQ           mixratio=Qa/(1-Qa)
# Line 185  cQQ           flwupa=-0.985*stefan*tsf** Line 184  cQQ           flwupa=-0.985*stefan*tsf**
184  cQQ  &                  *(0.39-0.05*sqrt(ea))  cQQ  &                  *(0.39-0.05*sqrt(ea))
185  cQQ  &                  *(1-0.6*nc**2)  cQQ  &                  *(1-0.6*nc**2)
186                if (iceornot.eq.0) then                if (iceornot.eq.0) then
187                 flwupa=-ocean_emissivity*stefan*tsf**4                 flwupa=ocean_emissivity*stefan*tsf**4
188                   dflwupdT=4. _d 0*ocean_emissivity*stefan*tsf**3
189                else                else
190                 if (iceornot.eq.2) then                 if (iceornot.eq.2) then
191                   flwupa=-snow_emissivity*stefan*tsf**4                  flwupa=snow_emissivity*stefan*tsf**4
192                    dflwupdT=4. _d 0*snow_emissivity*stefan*tsf**3
193                 else                 else
194                   flwupa=-ice_emissivity*stefan*tsf**4                  flwupa=ice_emissivity*stefan*tsf**4
195                    dflwupdT=4. _d 0*ice_emissivity*stefan*tsf**3
196                 endif                 endif
197                endif                endif
 #ifdef ALLOW_THERM_SEAICE  
 cswdcou -- remove ---  
 cQQ    if (iceornot.gt.0) then  
 cswdcou --- end remove ---  
 c derivatives with respect to Tsf  
198  cQQ           dflhdT = -clha*Tf0kel*ssq*22.47/(tsf**2)  cQQ           dflhdT = -clha*Tf0kel*ssq*22.47/(tsf**2)
199  c             dflhdT = -clha*Lath*ssq/(Rvap*tsf**2)  c             dflhdT = -clha*Lath*ssq/(Rvap*tsf**2)
200                dflhdT = -clha*ssq*Lath*2.166847E-3/(Tsf**2)                dflhdT = -clha*ssq*Lath*2.166847 _d -3/(Tsf**2)
201                dfshdT = -csha                dfshdT = -csha
202                if (iceornot.eq.0) then  cQQ           dflwupdT= 4.*0.985*stefan*tsf**3
                dflwupdT=-4.*ocean_emissivity*stefan*tsf**3  
               else  
                if (iceornot.eq.2) then  
                 dflwupdT=-4.*snow_emissivity*stefan*tsf**3  
                else  
                 dflwupdT=-4.*ice_emissivity*stefan*tsf**3  
                endif  
               endif  
 cQQ           dflwupdT=-4.*0.985*stefan*tsf**3  
203  cQQ  &                  *(0.39-0.05*sqrt(ea))  cQQ  &                  *(0.39-0.05*sqrt(ea))
204  cQQ  &                  *(1-0.6*nc**2)  cQQ  &                  *(1-0.6*nc**2)
205  c total derivative with respect to surface temperature  c total derivative with respect to surface temperature
206                df0dT=dflwupdT+dfshdT+dflhdT                df0dT=-dflwupdT+dfshdT+dflhdT
 cBB  
 cBB     print*,'derivatives:',df0dT,dflwupdT,dfshdT,dflhdT  
 cswdcou -- remove ---  
 cQQ     endif  
 cswdcou  
 #endif /*ALLOW_THERM_SEAICE*/  
207  c  c
208  c wind stress at center points  c wind stress at center points
209                if (.NOT.windread) then                if (.NOT.windread) then
# Line 230  c wind stress at center points Line 212  c wind stress at center points
212                 endif                 endif
213  #endif /*ALLOW_BULK_FORCE*/  #endif /*ALLOW_BULK_FORCE*/
214    
215        return        RETURN
216        end        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22