pkg/bulk_force/bulkf_formula_lanl.F

C $Header: /u/gcmpack/MITgcm/pkg/bulk_force/bulkf_formula_lanl.F,v 1.6 2005/11/04 01:28:39 jmc Exp $
C $Name:  $

#include "BULK_FORCE_OPTIONS.h"

      subroutine bulkf_formula_lanl(
     I                           uw, vw, us, Ta, Qa, nc, tsf_in,
     O                           flwupa, flha, fsha, df0dT,
     O                           ust, vst, evp, ssq, dEvdT, 
     I                           iceornot, myThid )

c swd -- bulkf formula used in bulkf and ice pkgs
c taken from exf package 

       IMPLICIT NONE
c
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "BULKF_PARAMS.h"

c
c calculate bulk formula fluxes over open ocean or seaice
c
c input
      _RL us                 ! wind speed
      _RL uw                 ! zonal wind speed (at grid center)
      _RL vw                 ! meridional wind speed (at grid center)
      _RL Ta                 ! air temperature at ht
      _RL Qa                 ! specific humidity at ht
      _RL tsf_in             ! surface temperature (either ice or open water)
      _RL nc                 ! fraction cloud cover
      integer iceornot       ! 0=open water, 1=ice cover
      integer myThid         ! my Thread Id number
c output
      _RL flwupa             ! upward long wave radiation (>0 upward)
      _RL flha               ! latent heat flux         (>0 downward)
      _RL fsha               ! sensible heat flux       (>0 downward)
      _RL df0dT              ! derivative of heat flux with respect to Tsf
      _RL ust                ! zonal wind stress (at grid center)
      _RL vst                ! meridional wind stress (at grid center)
      _RL evp                ! evaporation rate (over open water) [kg/m2/s]
      _RL ssq                ! surface specific humidity (kg/kg)
      _RL dEvdT              ! derivative of evap. with respect to Tsf [kg/m2/s/K]

#ifdef ALLOW_BULK_FORCE

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
      _RL dflhdT             ! derivative of latent heat with respect to T
      _RL dfshdT             ! derivative of sensible heat with respect to T
      _RL dflwupdT           ! derivative of long wave with respect to T
c     _RL mixratio
c     _RL ea
      _RL psim_fac
      _RL umin
      _RL lath
      _RL csha
      _RL clha
      _RL zice
      _RL ssq0, ssq1, ssq2   ! constant used in surface specific humidity
      _RL bulkf_Cdn          ! drag coefficient
      integer niter_bulk, iter

c --- external functions ---
c     _RL       exf_BulkCdn
c     external  exf_BulkCdn
c     _RL       exf_BulkqSat
c     external  exf_BulkqSat
c     _RL       exf_BulkRhn
c     external  exf_BulkRhn

      DATA   ssq0,           ssq1,           ssq2 
     &     / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /

cQQQQQQQQQQQQQQQQQQQQq
c -- compute turbulent surface fluxes
              ht =  2. _d 0
              hq =  2. _d 0
              hu = 10. _d 0
              zref = 10. _d 0
              zice = 0.0005 _d 0
              aln = log(ht/zref)
              niter_bulk = 5
              cdalton = 0.0346000 _d 0
              czol = zref*xkar*gravity
              psim_fac=5. _d 0
              umin=1. _d 0
c 
              lath=Lvap
              if (iceornot.gt.0) lath=Lvap+Lfresh
              Tsf=Tsf_in+Tf0kel
c     Wind speed 
              if (us.eq.0. _d 0) then
                us = sqrt(uw*uw + vw*vw)
              endif
              usm = max(us,umin)
cQQQ try to limit drag
cQQ           usm = min(usm,5. _d 0)
c
              t0     = Ta*(1. _d 0 + humid_fac*Qa)
cQQ           ssq    = 0.622*6.11*exp(22.47*(1.d0-Tf0kel/tsf))/p0
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
c             ssq = 3.797915 _d 0*exp(
c    &                lath*(7.93252 _d -6 - 2.166847 _d -3/Tsf)
c    &                               )/p0
              ssq = ssq0*exp( lath*(ssq1-ssq2/Tsf) ) / p0
cBB debugging
cBB             print*,'ice, ssq',  ssq
c
              deltap = ta  - tsf + gamma_blk*ht
              delq   = Qa - ssq
c
c initialize estimate exchange coefficients
              rdn=xkar/(log(zref/zice))
              rhn=rdn
              ren=rdn
c calculate turbulent scales
              ustar=rdn*usm
              tstar=rhn*deltap
              qstar=ren*delq
c
c interation with psi-functions to find transfer coefficients
              do iter=1,niter_bulk
                 huol   = czol/ustar**2 *(tstar/t0 +
     &                    qstar/(1. _d 0/humid_fac+Qa))
                 huol   = sign( min(abs(huol),10. _d 0), huol)
                 stable = 5. _d -1 + sign(5. _d -1 , huol)
                 xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*huol)),1. _d 0)
                 x      = sqrt(xsq)
                 psimh = -5. _d 0*huol*stable + (1. _d 0-stable)*
     &                    (2. _d 0*log(5. _d -1*(1. _d 0+x)) +
     &                     2. _d 0*log(5. _d -1*(1. _d 0+xsq)) -
     &                     2. _d 0*atan(x) + pi*.5 _d 0)
                 psixh  = -5. _d 0*huol*stable + (1. _d 0-stable)*
     &                     (2. _d 0*log(5. _d -1*(1. _d 0+xsq)))
   
c Update the transfer coefficients 
   
                 rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar)
                 rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar)
                 re = rh
c  Update ustar, tstar, qstar using updated, shifted coefficients.
                 ustar = rd*usm
                 qstar = re*delq
                 tstar = rh*deltap
              enddo
c
                tau   = rhoa*ustar**2
                tau   = tau*us/usm
                csha   = rhoa*cpair*us*rh*rd   
                clha   = rhoa*lath*us*re*rd   
c
                fsha  = csha*deltap
                flha  = clha*delq
                evp   = -flha/lath
c
c up long wave radiation
cQQ           mixratio=Qa/(1-Qa)
cQQ           ea=p0*mixratio/(0.62197+mixratio)
cQQ           flwupa=-0.985*stefan*tsf**4
cQQ  &                  *(0.39-0.05*sqrt(ea))
cQQ  &                  *(1-0.6*nc**2)
              if (iceornot.eq.0) then
               flwupa=ocean_emissivity*stefan*tsf**4
               dflwupdT=4. _d 0*ocean_emissivity*stefan*tsf**3
              else
               if (iceornot.eq.2) then
                flwupa=snow_emissivity*stefan*tsf**4
                dflwupdT=4. _d 0*snow_emissivity*stefan*tsf**3
               else 
                flwupa=ice_emissivity*stefan*tsf**4
                dflwupdT=4. _d 0*ice_emissivity*stefan*tsf**3
               endif
              endif
cQQ           dflhdT = -clha*Tf0kel*ssq*22.47/(tsf**2)
c             dflhdT = -clha*Lath*ssq/(Rvap*tsf**2)
c             dflhdT = -clha*ssq*Lath*2.166847 _d -3/(Tsf**2)
              dEvdT  =  clha*ssq*ssq2/(Tsf*Tsf)
              dflhdT = -lath*dEvdT
              dfshdT = -csha
cQQ           dflwupdT= 4.*0.985*stefan*tsf**3
cQQ  &                  *(0.39-0.05*sqrt(ea))
cQQ  &                  *(1-0.6*nc**2)
c total derivative with respect to surface temperature
              df0dT=-dflwupdT+dfshdT+dflhdT
c
c wind stress at center points
c             if (.NOT.windread) then
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
              bulkf_Cdn = cdrag_1/usm + cdrag_2 + cdrag_3*usm
              ust = rhoa*bulkf_Cdn*us*uw
              vst = rhoa*bulkf_Cdn*us*vw
#endif /*ALLOW_BULK_FORCE*/

      RETURN
      END
1	jmc	1.7	C $Header: /u/gcmpack/MITgcm/pkg/bulk_force/bulkf_formula_lanl.F,v 1.6 2005/11/04 01:28:39 jmc Exp $
2	edhill	1.3	C $Name: $
3	cheisey	1.1
4	edhill	1.3	#include "BULK_FORCE_OPTIONS.h"
5	jmc	1.4
6	cheisey	1.1	subroutine bulkf_formula_lanl(
7			I uw, vw, us, Ta, Qa, nc, tsf_in,
8	jmc	1.4	O flwupa, flha, fsha, df0dT,
9	jmc	1.5	O ust, vst, evp, ssq, dEvdT,
10	jmc	1.7	I iceornot, myThid )
11	jmc	1.4
12			c swd -- bulkf formula used in bulkf and ice pkgs
13			c taken from exf package
14	cheisey	1.1
15			IMPLICIT NONE
16			c
17			#include "EEPARAMS.h"
18			#include "SIZE.h"
19			#include "PARAMS.h"
20	jmc	1.4	#include "BULKF_PARAMS.h"
21	cheisey	1.1
22			c
23			c calculate bulk formula fluxes over open ocean or seaice
24			c
25			c input
26			_RL us ! wind speed
27			_RL uw ! zonal wind speed (at grid center)
28			_RL vw ! meridional wind speed (at grid center)
29			_RL Ta ! air temperature at ht
30			_RL Qa ! specific humidity at ht
31			_RL tsf_in ! surface temperature (either ice or open water)
32			_RL nc ! fraction cloud cover
33			integer iceornot ! 0=open water, 1=ice cover
34	jmc	1.7	integer myThid ! my Thread Id number
35	cheisey	1.1	c output
36	jmc	1.4	_RL flwupa ! upward long wave radiation (>0 upward)
37			_RL flha ! latent heat flux (>0 downward)
38			_RL fsha ! sensible heat flux (>0 downward)
39	cheisey	1.1	_RL df0dT ! derivative of heat flux with respect to Tsf
40			_RL ust ! zonal wind stress (at grid center)
41			_RL vst ! meridional wind stress (at grid center)
42	jmc	1.5	_RL evp ! evaporation rate (over open water) [kg/m2/s]
43	cheisey	1.1	_RL ssq ! surface specific humidity (kg/kg)
44	jmc	1.5	_RL dEvdT ! derivative of evap. with respect to Tsf [kg/m2/s/K]
45	jmc	1.4
46			#ifdef ALLOW_BULK_FORCE
47
48	cheisey	1.1	c local variables
49			_RL tsf ! surface temperature in K
50			_RL ht ! reference height (m)
51			_RL hq ! reference height for humidity (m)
52			_RL hu ! reference height for wind speed (m)
53			_RL zref ! reference height
54			_RL czol
55			_RL usm ! wind speed limited
56			_RL t0 ! virtual temperature (K)
57			_RL deltap ! potential temperature diff (K)
58			_RL delq ! specific humidity difference
59			_RL stable
60			_RL rdn ,ren, rhn
61			_RL ustar
62			_RL tstar
63			_RL qstar
64			_RL huol
65			_RL xsq
66			_RL x
67			_RL re
68			_RL rh
69			_RL tau
70			_RL psimh
71			_RL psixh
72			_RL rd
73			_RL aln
74			_RL cdalton
75			_RL dflhdT ! derivative of latent heat with respect to T
76			_RL dfshdT ! derivative of sensible heat with respect to T
77			_RL dflwupdT ! derivative of long wave with respect to T
78	jmc	1.6	c _RL mixratio
79			c _RL ea
80	cheisey	1.1	_RL psim_fac
81			_RL umin
82			_RL lath
83			_RL csha
84			_RL clha
85			_RL zice
86	jmc	1.5	_RL ssq0, ssq1, ssq2 ! constant used in surface specific humidity
87	jmc	1.7	_RL bulkf_Cdn ! drag coefficient
88	cheisey	1.1	integer niter_bulk, iter
89
90			c --- external functions ---
91	jmc	1.7	c _RL exf_BulkCdn
92			c external exf_BulkCdn
93	jmc	1.4	c _RL exf_BulkqSat
94			c external exf_BulkqSat
95			c _RL exf_BulkRhn
96			c external exf_BulkRhn
97	cheisey	1.1
98	jmc	1.5	DATA ssq0, ssq1, ssq2
99			& / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /
100
101	cheisey	1.1	cQQQQQQQQQQQQQQQQQQQQq
102			c -- compute turbulent surface fluxes
103	jmc	1.4	ht = 2. _d 0
104			hq = 2. _d 0
105			hu = 10. _d 0
106			zref = 10. _d 0
107			zice = 0.0005 _d 0
108	cheisey	1.1	aln = log(ht/zref)
109			niter_bulk = 5
110	jmc	1.4	cdalton = 0.0346000 _d 0
111	cheisey	1.1	czol = zrefxkargravity
112	jmc	1.4	psim_fac=5. _d 0
113			umin=1. _d 0
114	cheisey	1.1	c
115			lath=Lvap
116			if (iceornot.gt.0) lath=Lvap+Lfresh
117			Tsf=Tsf_in+Tf0kel
118			c Wind speed
119	jmc	1.4	if (us.eq.0. _d 0) then
120	cheisey	1.1	us = sqrt(uwuw + vwvw)
121			endif
122			usm = max(us,umin)
123			cQQQ try to limit drag
124	jmc	1.4	cQQ usm = min(usm,5. _d 0)
125	cheisey	1.1	c
126	jmc	1.4	t0 = Ta(1. _d 0 + humid_facQa)
127	cheisey	1.1	cQQ ssq = 0.6226.11exp(22.47*(1.d0-Tf0kel/tsf))/p0
128	jmc	1.4	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
129	jmc	1.5	c ssq = 3.797915 _d 0*exp(
130			c & lath*(7.93252 _d -6 - 2.166847 _d -3/Tsf)
131			c & )/p0
132			ssq = ssq0exp( lath(ssq1-ssq2/Tsf) ) / p0
133	cheisey	1.1	cBB debugging
134			cBB print*,'ice, ssq', ssq
135			c
136			deltap = ta - tsf + gamma_blk*ht
137			delq = Qa - ssq
138			c
139			c initialize estimate exchange coefficients
140			rdn=xkar/(log(zref/zice))
141			rhn=rdn
142			ren=rdn
143			c calculate turbulent scales
144			ustar=rdn*usm
145			tstar=rhn*deltap
146			qstar=ren*delq
147			c
148			c interation with psi-functions to find transfer coefficients
149			do iter=1,niter_bulk
150			huol = czol/ustar*2 (tstar/t0 +
151	jmc	1.4	& qstar/(1. _d 0/humid_fac+Qa))
152			huol = sign( min(abs(huol),10. _d 0), huol)
153	cheisey	1.1	stable = 5. _d -1 + sign(5. _d -1 , huol)
154	jmc	1.4	xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*huol)),1. _d 0)
155	cheisey	1.1	x = sqrt(xsq)
156	jmc	1.4	psimh = -5. _d 0huolstable + (1. _d 0-stable)*
157			& (2. _d 0log(5. _d -1(1. _d 0+x)) +
158			& 2. _d 0log(5. _d -1(1. _d 0+xsq)) -
159			& 2. _d 0atan(x) + pi.5 _d 0)
160			psixh = -5. _d 0huolstable + (1. _d 0-stable)*
161			& (2. _d 0log(5. _d -1(1. _d 0+xsq)))
162	cheisey	1.1
163			c Update the transfer coefficients
164
165	jmc	1.4	rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar)
166			rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar)
167	cheisey	1.1	re = rh
168			c Update ustar, tstar, qstar using updated, shifted coefficients.
169			ustar = rd*usm
170			qstar = re*delq
171			tstar = rh*deltap
172	jmc	1.4	enddo
173	cheisey	1.1	c
174			tau = rhoaustar*2
175			tau = tau*us/usm
176			csha = rhoacpairusrhrd
177			clha = rhoalathusrerd
178			c
179			fsha = csha*deltap
180			flha = clha*delq
181	jmc	1.5	evp = -flha/lath
182	cheisey	1.1	c
183			c up long wave radiation
184			cQQ mixratio=Qa/(1-Qa)
185			cQQ ea=p0*mixratio/(0.62197+mixratio)
186	cheisey	1.2	cQQ flwupa=-0.985stefantsf**4
187	cheisey	1.1	cQQ & (0.39-0.05sqrt(ea))
188			cQQ & (1-0.6nc**2)
189			if (iceornot.eq.0) then
190	jmc	1.4	flwupa=ocean_emissivitystefantsf**4
191			dflwupdT=4. _d 0ocean_emissivitystefantsf*3
192	cheisey	1.1	else
193			if (iceornot.eq.2) then
194	jmc	1.4	flwupa=snow_emissivitystefantsf**4
195			dflwupdT=4. _d 0snow_emissivitystefantsf*3
196	cheisey	1.1	else
197	jmc	1.4	flwupa=ice_emissivitystefantsf**4
198			dflwupdT=4. _d 0ice_emissivitystefantsf*3
199	cheisey	1.1	endif
200			endif
201			cQQ dflhdT = -clhaTf0kelssq22.47/(tsf*2)
202			c dflhdT = -clhaLathssq/(Rvaptsf*2)
203	jmc	1.5	c dflhdT = -clhassqLath2.166847 _d -3/(Tsf*2)
204			dEvdT = clhassqssq2/(Tsf*Tsf)
205			dflhdT = -lath*dEvdT
206	cheisey	1.1	dfshdT = -csha
207	jmc	1.4	cQQ dflwupdT= 4.0.985stefantsf*3
208	cheisey	1.1	cQQ & (0.39-0.05sqrt(ea))
209			cQQ & (1-0.6nc**2)
210			c total derivative with respect to surface temperature
211	jmc	1.4	df0dT=-dflwupdT+dfshdT+dflhdT
212	cheisey	1.1	c
213			c wind stress at center points
214	jmc	1.7	c if (.NOT.windread) then
215			c ust = rhoaexf_BulkCdn(usm)us*uw
216			c vst = rhoaexf_BulkCdn(usm)us*vw
217			c endif
218			C- In-lining of function: exf_BulkCdn(umps) = cdrag_1/umps + cdrag_2 + cdrag_3*umps
219			bulkf_Cdn = cdrag_1/usm + cdrag_2 + cdrag_3*usm
220			ust = rhoabulkf_Cdnus*uw
221			vst = rhoabulkf_Cdnus*vw
222	cheisey	1.1	#endif /ALLOW_BULK_FORCE/
223
224	jmc	1.4	RETURN
225			END