pkg/bulk_force/bulkf_formula_lanl.F

C $Header: /u/gcmpack/MITgcm/pkg/bulk_force/bulkf_formula_lanl.F,v 1.7 2006/01/22 15:51:35 jmc Exp $
C $Name:  $

#include "BULK_FORCE_OPTIONS.h"

CBOP
C     !ROUTINE: BULKF_FORMULA_LANL
C     !INTERFACE:
      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     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE BULKF_FORMULA_LANL
c     | o Calculate bulk formula fluxes over open ocean or seaice
C     *==========================================================*
C     \ev
C swd -- bulkf formula used in bulkf and ice pkgs
C        taken from exf package
C
C     wind stress = (ust,vst) = rhoA * Cd * Ws * (del.u,del.v)
C     Sensib Heat flux = fsha = rhoA * Ch * Ws * del.T * CpAir
C     Latent Heat flux = flha = rhoA * Ce * Ws * del.Q * Lvap
C                      = -Evap * Lvap
C   with Ws = wind speed = sqrt(del.u^2 +del.v^2) ;
C        del.T = Tair - Tsurf ; del.Q = Qair - Qsurf
C        Cd,Ch,Ce = transfer coefficient for momentum, sensible
C                    & latent heat flux [no units]
C     *==========================================================*

C     !USES:
       IMPLICIT NONE
C     === Global variables ===
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "BULKF_PARAMS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     input:
      _RL uw                 ! zonal wind speed (at grid center) [m/s]
      _RL vw                 ! meridional wind speed (at grid center) [m/s]
      _RL us                 ! wind speed        [m/s]   at height hu
      _RL Ta                 ! air temperature   [K]     at height ht
      _RL Qa                 ! specific humidity [kg/kg] at heigth ht
      _RL nc                 ! fraction cloud cover
      _RL tsf_in             ! sea-ice or sea surface temperature [oC]
      INTEGER iceornot       ! 0=open water, 1=sea-ice, 2=sea-ice with snow
      INTEGER myThid         ! my Thread Id number
C     output:
      _RL flwupa             ! upward long wave radiation (>0 upward) [W/m2]
      _RL flha               ! latent heat flux         (>0 downward) [W/m2]
      _RL fsha               ! sensible heat flux       (>0 downward) [W/m2]
      _RL df0dT              ! derivative of heat flux with respect to Tsf [W/m2/K]
      _RL ust                ! zonal wind stress (at grid center)     [N/m2]
      _RL vst                ! meridional wind stress (at grid center)[N/m2]
      _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]
CEOP

#ifdef ALLOW_BULK_FORCE

C     == Local variables ==
      _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 tsf                ! surface temperature [K]
      _RL ht                 ! height for air temperature [m]
c     _RL hq                 ! height for humidity [m]
      _RL hu                 ! height for wind speed [m]
c     _RL zref               ! reference height [m]
      _RL usm                ! wind speed limited [m/s]
c     _RL umin               ! minimum wind speed used for drag-coeff [m/s]
      _RL lath               ! latent heat of vaporization or sublimation
      _RL t0                 ! virtual temperature [K]
      _RL deltap             ! potential temperature diff [K]
      _RL delq               ! specific humidity difference [kg/kg]
      _RL ustar              ! friction velocity [m/s]
      _RL tstar              ! temperature scale [K]
      _RL qstar              ! humidity scale  [kg/kg]
      _RL rd                 ! = sqrt(Cd)          [-]
      _RL re                 ! = Ce / sqrt(Cd)     [-]
      _RL rh                 ! = Ch / sqrt(Cd)     [-]
      _RL rdn, ren, rhn      ! initial (neutral) values of rd, re, rh
      _RL stable             ! = 1 if stable ; = 0 if unstable
      _RL huol               ! stability parameter [-]
      _RL x                  ! stability function  [-]
      _RL xsq                ! = x^2               [-]
      _RL psimh              ! momentum stability function
      _RL psixh              ! latent & sensib. stability function
      _RL czol               ! = zref*Karman_cst*gravity
      _RL aln                ! = log(ht/zref)
c     _RL cdalton            ! coeff to evaluate Dalton Number
c     _RL mixratio
c     _RL ea
c     _RL psim_fac
      _RL tau                ! surface stress  coef = ?
      _RL csha               ! sensib.heat flx coef = rhoA * Ws * Ch * CpAir
      _RL clha               ! latent heat flx coef = rhoA * Ws * Ce * Lvap
      _RL zice
      _RL ssq0, ssq1, ssq2   ! constant used in surface specific humidity
      _RL p0                 ! reference sea-level atmospheric pressure [mb]
      _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 /
      DATA   p0 / 1013. _d 0 /

C--- Compute turbulent surface fluxes
              ht =  2. _d 0
c             hq =  2. _d 0
              hu = 10. _d 0
c             zref = 10. _d 0
              zice = 0.0005 _d 0
              aln = log(ht/zref)
              niter_bulk = 5
c             cdalton = 0.0346000 _d 0
              czol = zref*xkar*gravity
c             psim_fac=5. _d 0
c             umin=1. _d 0

              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)
c
              t0     = Ta*(1. _d 0 + humid_fac*Qa)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
cQQ           ssq    = 0.622*6.11*exp(22.47*(1.d0-Tf0kel/tsf))/p0
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

              deltap = ta  - tsf + gamma_blk*ht
              delq   = Qa - ssq

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

              tau   = rhoa*ustar**2
              tau   = tau*us/usm
              csha  = rhoa*cpair*us*rh*rd
              clha  = rhoa*lath*us*re*rd

              fsha  = csha*deltap
              flha  = clha*delq
              evp   = -flha/lath

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