pkg/bulk_force/bulkf_formula_lanl.F

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