pkg/bulk_force/bulkf_formula_lanl.F

C $Header:  $
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, 
     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)
      _RL ssq                ! surface specific humidity (kg/kg)

#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
      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

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