pkg/bulk_force/bulkf_formula_lanl.F

C $Header: /u/gcmpack/MITgcm/pkg/bulk_force/bulkf_formula_lanl.F,v 1.4 2003/11/23 01:36:55 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, windread )

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
      logical windread       !
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 htol
      _RL hqol
      _RL xsq
      _RL x
      _RL re
      _RL rh
      _RL tau
      _RL psimh
      _RL psixh
      _RL rd
      _RL uzn
      _RL shn
      _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
      _RL mixratio
      _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
      integer niter_bulk, iter

c --- external functions ---
      _RL       exf_BulkCdn
      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
              if (.NOT.windread) then
                 ust = rhoa*exf_BulkCdn(usm)*us*uw
                 vst = rhoa*exf_BulkCdn(usm)*us*vw
               endif
#endif /*ALLOW_BULK_FORCE*/

      RETURN
      END
1	jmc	1.5	C $Header: /u/gcmpack/MITgcm/pkg/bulk_force/bulkf_formula_lanl.F,v 1.4 2003/11/23 01:36:55 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.4	I iceornot, windread )
11
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			logical windread !
35			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 htol
66			_RL hqol
67			_RL xsq
68			_RL x
69			_RL re
70			_RL rh
71			_RL tau
72			_RL psimh
73			_RL psixh
74			_RL rd
75			_RL uzn
76			_RL shn
77			_RL aln
78			_RL cdalton
79			_RL dflhdT ! derivative of latent heat with respect to T
80			_RL dfshdT ! derivative of sensible heat with respect to T
81			_RL dflwupdT ! derivative of long wave with respect to T
82			_RL mixratio
83			_RL ea
84			_RL psim_fac
85			_RL umin
86			_RL lath
87			_RL csha
88			_RL clha
89			_RL zice
90	jmc	1.5	_RL ssq0, ssq1, ssq2 ! constant used in surface specific humidity
91	cheisey	1.1	integer niter_bulk, iter
92
93			c --- external functions ---
94			_RL exf_BulkCdn
95			external exf_BulkCdn
96	jmc	1.4	c _RL exf_BulkqSat
97			c external exf_BulkqSat
98			c _RL exf_BulkRhn
99			c external exf_BulkRhn
100	cheisey	1.1
101	jmc	1.5	DATA ssq0, ssq1, ssq2
102			& / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /
103
104	cheisey	1.1	cQQQQQQQQQQQQQQQQQQQQq
105			c -- compute turbulent surface fluxes
106	jmc	1.4	ht = 2. _d 0
107			hq = 2. _d 0
108			hu = 10. _d 0
109			zref = 10. _d 0
110			zice = 0.0005 _d 0
111	cheisey	1.1	aln = log(ht/zref)
112			niter_bulk = 5
113	jmc	1.4	cdalton = 0.0346000 _d 0
114	cheisey	1.1	czol = zrefxkargravity
115	jmc	1.4	psim_fac=5. _d 0
116			umin=1. _d 0
117	cheisey	1.1	c
118			lath=Lvap
119			if (iceornot.gt.0) lath=Lvap+Lfresh
120			Tsf=Tsf_in+Tf0kel
121			c Wind speed
122	jmc	1.4	if (us.eq.0. _d 0) then
123	cheisey	1.1	us = sqrt(uwuw + vwvw)
124			endif
125			usm = max(us,umin)
126			cQQQ try to limit drag
127	jmc	1.4	cQQ usm = min(usm,5. _d 0)
128	cheisey	1.1	c
129	jmc	1.4	t0 = Ta(1. _d 0 + humid_facQa)
130	cheisey	1.1	cQQ ssq = 0.6226.11exp(22.47*(1.d0-Tf0kel/tsf))/p0
131	jmc	1.4	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
132	jmc	1.5	c ssq = 3.797915 _d 0*exp(
133			c & lath*(7.93252 _d -6 - 2.166847 _d -3/Tsf)
134			c & )/p0
135			ssq = ssq0exp( lath(ssq1-ssq2/Tsf) ) / p0
136	cheisey	1.1	cBB debugging
137			cBB print*,'ice, ssq', ssq
138			c
139			deltap = ta - tsf + gamma_blk*ht
140			delq = Qa - ssq
141			c
142			c initialize estimate exchange coefficients
143			rdn=xkar/(log(zref/zice))
144			rhn=rdn
145			ren=rdn
146			c calculate turbulent scales
147			ustar=rdn*usm
148			tstar=rhn*deltap
149			qstar=ren*delq
150			c
151			c interation with psi-functions to find transfer coefficients
152			do iter=1,niter_bulk
153			huol = czol/ustar*2 (tstar/t0 +
154	jmc	1.4	& qstar/(1. _d 0/humid_fac+Qa))
155			huol = sign( min(abs(huol),10. _d 0), huol)
156	cheisey	1.1	stable = 5. _d -1 + sign(5. _d -1 , huol)
157	jmc	1.4	xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*huol)),1. _d 0)
158	cheisey	1.1	x = sqrt(xsq)
159	jmc	1.4	psimh = -5. _d 0huolstable + (1. _d 0-stable)*
160			& (2. _d 0log(5. _d -1(1. _d 0+x)) +
161			& 2. _d 0log(5. _d -1(1. _d 0+xsq)) -
162			& 2. _d 0atan(x) + pi.5 _d 0)
163			psixh = -5. _d 0huolstable + (1. _d 0-stable)*
164			& (2. _d 0log(5. _d -1(1. _d 0+xsq)))
165	cheisey	1.1
166			c Update the transfer coefficients
167
168	jmc	1.4	rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar)
169			rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar)
170	cheisey	1.1	re = rh
171			c Update ustar, tstar, qstar using updated, shifted coefficients.
172			ustar = rd*usm
173			qstar = re*delq
174			tstar = rh*deltap
175	jmc	1.4	enddo
176	cheisey	1.1	c
177			tau = rhoaustar*2
178			tau = tau*us/usm
179			csha = rhoacpairusrhrd
180			clha = rhoalathusrerd
181			c
182			fsha = csha*deltap
183			flha = clha*delq
184	jmc	1.5	evp = -flha/lath
185	cheisey	1.1	c
186			c up long wave radiation
187			cQQ mixratio=Qa/(1-Qa)
188			cQQ ea=p0*mixratio/(0.62197+mixratio)
189	cheisey	1.2	cQQ flwupa=-0.985stefantsf**4
190	cheisey	1.1	cQQ & (0.39-0.05sqrt(ea))
191			cQQ & (1-0.6nc**2)
192			if (iceornot.eq.0) then
193	jmc	1.4	flwupa=ocean_emissivitystefantsf**4
194			dflwupdT=4. _d 0ocean_emissivitystefantsf*3
195	cheisey	1.1	else
196			if (iceornot.eq.2) then
197	jmc	1.4	flwupa=snow_emissivitystefantsf**4
198			dflwupdT=4. _d 0snow_emissivitystefantsf*3
199	cheisey	1.1	else
200	jmc	1.4	flwupa=ice_emissivitystefantsf**4
201			dflwupdT=4. _d 0ice_emissivitystefantsf*3
202	cheisey	1.1	endif
203			endif
204			cQQ dflhdT = -clhaTf0kelssq22.47/(tsf*2)
205			c dflhdT = -clhaLathssq/(Rvaptsf*2)
206	jmc	1.5	c dflhdT = -clhassqLath2.166847 _d -3/(Tsf*2)
207			dEvdT = clhassqssq2/(Tsf*Tsf)
208			dflhdT = -lath*dEvdT
209	cheisey	1.1	dfshdT = -csha
210	jmc	1.4	cQQ dflwupdT= 4.0.985stefantsf*3
211	cheisey	1.1	cQQ & (0.39-0.05sqrt(ea))
212			cQQ & (1-0.6nc**2)
213			c total derivative with respect to surface temperature
214	jmc	1.4	df0dT=-dflwupdT+dfshdT+dflhdT
215	cheisey	1.1	c
216			c wind stress at center points
217			if (.NOT.windread) then
218			ust = rhoaexf_BulkCdn(usm)us*uw
219			vst = rhoaexf_BulkCdn(usm)us*vw
220			endif
221			#endif /ALLOW_BULK_FORCE/
222
223	jmc	1.4	RETURN
224			END