New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sbcblk.F90 in NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90 @ 11275

Last change on this file since 11275 was 11275, checked in by smasson, 5 years ago
dev_r11265_ABL : add ABL interface, see #2131
Property svn:keywords set to `Id`
File size: 74.7 KB

Line
1	MODULE sbcblk
2	!!======================================================================
3	!! * MODULE sbcblk *
4	!! Ocean forcing: momentum, heat and freshwater flux formulation
5	!! Aerodynamic Bulk Formulas
6	!! SUCCESSOR OF "sbcblk_core"
7	!!=====================================================================
8	!! History : 1.0 ! 2004-08 (U. Schweckendiek) Original CORE code
9	!! 2.0 ! 2005-04 (L. Brodeau, A.M. Treguier) improved CORE bulk and its user interface
10	!! 3.0 ! 2006-06 (G. Madec) sbc rewritting
11	!! - ! 2006-12 (L. Brodeau) Original code for turb_core
12	!! 3.2 ! 2009-04 (B. Lemaire) Introduce iom_put
13	!! 3.3 ! 2010-10 (S. Masson) add diurnal cycle
14	!! 3.4 ! 2011-11 (C. Harris) Fill arrays required by CICE
15	!! 3.7 ! 2014-06 (L. Brodeau) simplification and optimization of CORE bulk
16	!! 4.0 ! 2016-06 (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
17	!! ! ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
18	!! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine
19	!! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle)
20	!! 4.0 ! 2019-03 (F. Lemarié & G. Samson) add ABL compatibility (ln_abl=TRUE)
21	!!----------------------------------------------------------------------
22
23	!!----------------------------------------------------------------------
24	!! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition
25	!! sbc_blk : bulk formulation as ocean surface boundary condition
26	!! blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model (ln_abl=TRUE)
27	!! blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step (ln_abl=TRUE)
28	!! rho_air : density of (moist) air (depends on T_air, q_air and SLP
29	!! cp_air : specific heat of (moist) air (depends spec. hum. q_air)
30	!! q_sat : saturation humidity as a function of SLP and temperature
31	!! L_vap : latent heat of vaporization of water as a function of temperature
32	!! sea-ice case only :
33	!! blk_ice_tau : provide the air-ice stress
34	!! blk_ice_flx : provide the heat and mass fluxes at air-ice interface
35	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
36	!! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
37	!! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
38	!!----------------------------------------------------------------------
39	USE oce ! ocean dynamics and tracers
40	USE dom_oce ! ocean space and time domain
41	USE phycst ! physical constants
42	USE fldread ! read input fields
43	USE sbc_oce ! Surface boundary condition: ocean fields
44	USE cyclone ! Cyclone 10m wind form trac of cyclone centres
45	USE sbcdcy ! surface boundary condition: diurnal cycle
46	USE sbcwave , ONLY : cdn_wave ! wave module
47	USE sbc_ice ! Surface boundary condition: ice fields
48	USE lib_fortran ! to use key_nosignedzero
49	#if defined key_si3
50	USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif
51	USE icethd_dh ! for CALL ice_thd_snwblow
52	#endif
53	USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009)
54	USE sbcblk_algo_coare ! => turb_coare : COAREv3.0 (Fairall et al. 2003)
55	USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)
56	USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 31)
57	!
58	USE iom ! I/O manager library
59	USE in_out_manager ! I/O manager
60	USE lib_mpp ! distribued memory computing library
61	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
62	USE prtctl ! Print control
63
64	IMPLICIT NONE
65	PRIVATE
66
67	PUBLIC sbc_blk_init ! called in sbcmod
68	PUBLIC sbc_blk ! called in sbcmod
69	PUBLIC blk_oce_1 ! called in sbcabl
70	PUBLIC blk_oce_2 ! called in sbcabl
71	PUBLIC rho_air ! called in ablmod
72	PUBLIC cp_air ! called in ablmod
73	#if defined key_si3
74	PUBLIC blk_ice_tau ! routine called in icesbc
75	PUBLIC blk_ice_flx ! routine called in icesbc
76	PUBLIC blk_ice_qcn ! routine called in icesbc
77	#endif
78
79	INTERFACE cp_air
80	MODULE PROCEDURE cp_air_0d, cp_air_2d
81	END INTERFACE
82
83	!!Lolo: should ultimately be moved in the module with all physical constants ?
84	!!gm : In principle, yes.
85	REAL(wp), PARAMETER :: Cp_dry = 1005.0 !: Specic heat of dry air, constant pressure [J/K/kg]
86	REAL(wp), PARAMETER :: Cp_vap = 1860.0 !: Specic heat of water vapor, constant pressure [J/K/kg]
87	REAL(wp), PARAMETER :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg]
88	REAL(wp), PARAMETER :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg]
89	REAL(wp), PARAMETER :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
90	REAL(wp), PARAMETER :: rctv0 = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
91
92	INTEGER , PUBLIC, PARAMETER :: jpfld =11 !: maximum number of files to read
93	INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 !: index of 10m wind velocity (i-component) (m/s) at T-point
94	INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 !: index of 10m wind velocity (j-component) (m/s) at T-point
95	INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 !: index of 10m air temperature (Kelvin)
96	INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 !: index of specific humidity ( % )
97	INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 !: index of solar heat (W/m2)
98	INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 !: index of Long wave (W/m2)
99	INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 !: index of total precipitation (rain+snow) (Kg/m2/s)
100	INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 !: index of snow (solid prcipitation) (kg/m2/s)
101	INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 !: index of sea level pressure (Pa)
102	INTEGER , PUBLIC, PARAMETER :: jp_hpgi =10 !: index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point
103	INTEGER , PUBLIC, PARAMETER :: jp_hpgj =11 !: index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point
104
105	TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
106
107	! !!! Bulk parameters
108	REAL(wp), PARAMETER :: cpa = 1000.5 ! specific heat of air (only used for ice fluxes now...)
109	REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation
110	REAL(wp), PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant
111	REAL(wp), PARAMETER :: Cd_ice = 1.4e-3 ! transfer coefficient over ice
112	REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant
113	!
114	! !!* Namelist namsbc_blk : bulk parameters
115	LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008)
116	LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003)
117	LOGICAL :: ln_COARE_3p5 ! "COARE 3.5" algorithm (Edson et al. 2013)
118	LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 31)
119	!
120	REAL(wp) :: rn_pfac ! multiplication factor for precipitation
121	REAL(wp) :: rn_efac ! multiplication factor for evaporation
122	REAL(wp) :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress
123	REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements
124	REAL(wp) :: rn_zu ! z(u) : height of wind measurements
125	!!gm ref namelist initialize it so remove the setting to false below
126	LOGICAL :: ln_Cd_L12 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012)
127	LOGICAL :: ln_Cd_L15 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
128	!
129	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_atm ! transfer coefficient for momentum (tau)
130	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ch_atm ! transfer coefficient for sensible heat (Q_sens)
131	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ce_atm ! tansfert coefficient for evaporation (Q_lat)
132	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_zu ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme)
133	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_zu ! air spec. hum. at wind speed height (needed by Lupkes 2015 bulk scheme)
134	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme
135
136	INTEGER :: nblk ! choice of the bulk algorithm
137	! ! associated indices:
138	INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008)
139	INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003)
140	INTEGER, PARAMETER :: np_COARE_3p5 = 3 ! "COARE 3.5" algorithm (Edson et al. 2013)
141	INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 31)
142
143	!! * Substitutions
144	# include "vectopt_loop_substitute.h90"
145	!!----------------------------------------------------------------------
146	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
147	!! $Id$
148	!! Software governed by the CeCILL license (see ./LICENSE)
149	!!----------------------------------------------------------------------
150	CONTAINS
151
152	INTEGER FUNCTION sbc_blk_alloc()
153	!!-------------------------------------------------------------------
154	!! * ROUTINE sbc_blk_alloc *
155	!!-------------------------------------------------------------------
156	ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
157	& cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
158	!
159	CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
160	IF( sbc_blk_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' )
161	END FUNCTION sbc_blk_alloc
162
163
164	SUBROUTINE sbc_blk_init
165	!!---------------------------------------------------------------------
166	!! * ROUTINE sbc_blk_init *
167	!!
168	!! ** Purpose : choose and initialize a bulk formulae formulation
169	!!
170	!! ** Method :
171	!!
172	!!----------------------------------------------------------------------
173	INTEGER :: ifpr ! dummy loop indice and argument
174	INTEGER :: ios, ierror, ioptio ! Local integer
175	!!
176	CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files
177	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
178	TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read
179	TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " "
180	TYPE(FLD_N) :: sn_slp , sn_hpgi, sn_hpgj ! " "
181	NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields
182	& sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj, &
183	& ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm
184	& cn_dir , rn_zqt, rn_zu, &
185	& rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15
186	!!---------------------------------------------------------------------
187	!
188	! ! allocate sbc_blk_core array
189	IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
190	!
191	! !** read bulk namelist
192	REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters
193	READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
194	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist', lwp )
195	!
196	REWIND( numnam_cfg ) !* Namelist namsbc_blk in configuration namelist : bulk parameters
197	READ ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
198	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist', lwp )
199	!
200	IF(lwm) WRITE( numond, namsbc_blk )
201	!
202	! !** initialization of the chosen bulk formulae (+ check)
203	! !* select the bulk chosen in the namelist and check the choice
204	ioptio = 0
205	IF( ln_NCAR ) THEN ; nblk = np_NCAR ; ioptio = ioptio + 1 ; ENDIF
206	IF( ln_COARE_3p0 ) THEN ; nblk = np_COARE_3p0 ; ioptio = ioptio + 1 ; ENDIF
207	IF( ln_COARE_3p5 ) THEN ; nblk = np_COARE_3p5 ; ioptio = ioptio + 1 ; ENDIF
208	IF( ln_ECMWF ) THEN ; nblk = np_ECMWF ; ioptio = ioptio + 1 ; ENDIF
209	!
210	IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
211	!
212	IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr
213	IF( sn_qsr%freqh /= 24. ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
214	IF( sn_qsr%ln_tint ) THEN
215	CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', &
216	& ' ==> We force time interpolation = .false. for qsr' )
217	sn_qsr%ln_tint = .false.
218	ENDIF
219	ENDIF
220	! !* set the bulk structure
221	! !- store namelist information in an array
222	slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj
223	slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw
224	slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi
225	slf_i(jp_prec) = sn_prec ; slf_i(jp_snow) = sn_snow
226	slf_i(jp_slp) = sn_slp
227	slf_i(jp_hpgi) = sn_hpgi ; slf_i(jp_hpgj) = sn_hpgj
228	!
229	! !- allocate the bulk structure
230	ALLOCATE( sf(jpfld), STAT=ierror )
231	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
232	! !- fill the bulk structure with namelist informations
233	CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
234	!
235	DO ifpr= 1, jpfld
236	!
237	IF( TRIM(sf(ifpr)%clrootname) /= 'NOT USED' ) THEN
238	IF( ln_abl .AND. ( ifpr == jp_wndi .OR. ifpr == jp_wndj .OR. ifpr == jp_humi &
239	& .OR. ifpr == jp_hpgi .OR. ifpr == jp_hpgj .OR. ifpr == jp_tair ) ) THEN
240	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,jpka) )
241	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,jpka,2) )
242	ELSE
243	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
244	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
245	ENDIF
246	!
247	IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 ) &
248	& CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', &
249	& ' This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
250	ENDIF
251	END DO
252	!
253	IF ( ln_wave ) THEN
254	!Activated wave module but neither drag nor stokes drift activated
255	IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN
256	CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
257	!drag coefficient read from wave model definable only with mfs bulk formulae and core
258	ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR ) THEN
259	CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
260	ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN
261	CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
262	ENDIF
263	ELSE
264	IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) &
265	& CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', &
266	& 'with drag coefficient (ln_cdgw =T) ' , &
267	& 'or Stokes Drift (ln_sdw=T) ' , &
268	& 'or ocean stress modification due to waves (ln_tauwoc=T) ', &
269	& 'or Stokes-Coriolis term (ln_stcori=T)' )
270	ENDIF
271	!
272	IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature and humidity
273	rn_zqt = ght_abl(2) ! set the bulk altitude to ABL first level
274	rn_zu = ght_abl(2)
275	IF(lwp) WRITE(numout,*)
276	IF(lwp) WRITE(numout,*) ' ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude'
277	ENDIF
278	!
279	IF(lwp) THEN !** Control print
280	!
281	WRITE(numout,) ! namelist
282	WRITE(numout,*) ' Namelist namsbc_blk (other than data information):'
283	WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR
284	WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0
285	WRITE(numout,*) ' "COARE 3.5" algorithm (Edson et al. 2013) ln_COARE_3p5 = ', ln_COARE_3p0
286	WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 31) ln_ECMWF = ', ln_ECMWF
287	WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt
288	WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu
289	WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac
290	WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac
291	WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac
292	WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))'
293	WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12
294	WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15
295	!
296	WRITE(numout,*)
297	SELECT CASE( nblk ) !* Print the choice of bulk algorithm
298	CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)'
299	CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)'
300	CASE( np_COARE_3p5 ) ; WRITE(numout,*) ' ==>>> "COARE 3.5" algorithm (Edson et al. 2013)'
301	CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 31)'
302	END SELECT
303	!
304	ENDIF
305	!
306	END SUBROUTINE sbc_blk_init
307
308
309	SUBROUTINE sbc_blk( kt )
310	!!---------------------------------------------------------------------
311	!! * ROUTINE sbc_blk *
312	!!
313	!! ** Purpose : provide at each time step the surface ocean fluxes
314	!! (momentum, heat, freshwater and runoff)
315	!!
316	!! ** Method : (1) READ each fluxes in NetCDF files:
317	!! the 10m wind velocity (i-component) (m/s) at T-point
318	!! the 10m wind velocity (j-component) (m/s) at T-point
319	!! the 10m or 2m specific humidity (kg/kg)
320	!! the solar heat (W/m2)
321	!! the Long wave (W/m2)
322	!! the 10m or 2m air temperature (Kelvin)
323	!! the total precipitation (rain+snow) (Kg/m2/s)
324	!! the snow (solid precipitation) (kg/m2/s)
325	!! ABL dynamical forcing (i/j-components of either hpg or geostrophic winds)
326	!! (2) CALL blk_oce_1 and blk_oce_2
327	!!
328	!! C A U T I O N : never mask the surface stress fields
329	!! the stress is assumed to be in the (i,j) mesh referential
330	!!
331	!! ** Action : defined at each time-step at the air-sea interface
332	!! - utau, vtau i- and j-component of the wind stress
333	!! - taum wind stress module at T-point
334	!! - wndm wind speed module at T-point over free ocean or leads in presence of sea-ice
335	!! - qns, qsr non-solar and solar heat fluxes
336	!! - emp upward mass flux (evapo. - precip.)
337	!! - sfx salt flux due to freezing/melting (non-zero only if ice is present)
338	!!
339	!! ** References : Large & Yeager, 2004 / Large & Yeager, 2008
340	!! Brodeau et al. Ocean Modelling 2010
341	!!----------------------------------------------------------------------
342	INTEGER, INTENT(in) :: kt ! ocean time step
343	!!
344	REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zevp
345	!!---------------------------------------------------------------------
346	!
347	CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step
348	!
349	! ! compute the surface ocean fluxes using bulk formulea
350	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
351	!
352	CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), & ! <<= in
353	& sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), & ! <<= in
354	& sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in
355	& zssq, zcd_du, zsen, zevp ) ! =>> out
356
357	CALL blk_oce_2( kt, sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1), & ! <<= in
358	& sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), & ! <<= in
359	& sf(jp_snow)%fnow(:,:,1), sst_m, & ! <<= in
360	& zsen, zevp ) ! <=> in out
361
362	ENDIF
363
364	#if defined key_cice
365	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
366	qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1)
367	IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
368	ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1)
369	ENDIF
370	tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)
371	qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1)
372	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
373	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
374	wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
375	wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
376	ENDIF
377	#endif
378	!
379	END SUBROUTINE sbc_blk
380
381
382	SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi, & ! inp
383	& pslp, pst , pu , pv, & ! inp
384	& pssq, pcd_du, psen, pevp ) ! out
385	!!---------------------------------------------------------------------
386	!! * ROUTINE blk_oce_1 *
387	!!
388	!! ** Purpose : Computes Cd x \|U\|, Ch x \|U\|, Ce x \|U\| for ABL integration
389	!!
390	!! ** Method : bulk formulae using atmospheric
391	!! fields from the ABL model at previous time-step
392	!!
393	!! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg)
394	!! - pcd_du : Cd x \|dU\| at T-points (m/s)
395	!! - psen : Ch x \|dU\| at T-points (m/s)
396	!! - pevp : Ce x \|dU\| at T-points (m/s)
397	!!---------------------------------------------------------------------
398	INTEGER , INTENT(in ) :: kt ! time step index
399	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at U-point [m/s]
400	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndj ! atmospheric wind at V-point [m/s]
401	REAL(wp), INTENT(in ), DIMENSION(:,:) :: phumi ! specific humidity at T-points [kg/kg]
402	REAL(wp), INTENT(in ), DIMENSION(:,:) :: ptair ! potential temperature at T-points [Kelvin]
403	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pslp ! sea-level pressure [Pa]
404	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pst ! surface temperature [Celcius]
405	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s]
406	REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s]
407	REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg]
408	REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du ! Cd x \|dU\| at T-points [m/s]
409	REAL(wp), INTENT( out), DIMENSION(:,:) :: psen ! Ch x \|dU\| at T-points [m/s]
410	REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp ! Ce x \|dU\| at T-points [m/s]
411	!
412	INTEGER :: ji, jj ! dummy loop indices
413	REAL(wp) :: zztmp ! local variable
414	REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point
415	REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin
416	REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]
417	REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K]
418	REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! density of air [kg/m^3]
419	!!---------------------------------------------------------------------
420	!
421	! local scalars ( place there for vector optimisation purposes)
422	zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
423
424	! ----------------------------------------------------------------------------- !
425	! 0 Wind components and module at T-point relative to the moving ocean !
426	! ----------------------------------------------------------------------------- !
427
428	! ... components ( U10m - U_oce ) at T-point (unmasked)
429	#if defined key_cyclone
430	zwnd_i(:,:) = 0._wp
431	zwnd_j(:,:) = 0._wp
432	CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012)
433	DO jj = 2, jpjm1
434	DO ji = fs_2, fs_jpim1 ! vect. opt.
435	pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
436	pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj)
437	END DO
438	END DO
439	#endif
440	DO jj = 2, jpjm1
441	DO ji = fs_2, fs_jpim1 ! vect. opt.
442	zwnd_i(ji,jj) = ( pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) )
443	zwnd_j(ji,jj) = ( pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) )
444	END DO
445	END DO
446	CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
447	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
448	wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) &
449	& + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1)
450
451	! ----------------------------------------------------------------------------- !
452	! I Turbulent FLUXES !
453	! ----------------------------------------------------------------------------- !
454
455	! ... specific humidity at SST and IST tmask(
456	pssq(:,:) = 0.98 * q_sat( zst(:,:), pslp(:,:) )
457	!
458	IF( ln_abl ) THEN
459	ztpot = ptair(:,:)
460	ELSE
461	! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
462	! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
463	! (since reanalysis products provide T at z, not theta !)
464	ztpot = ptair(:,:) + gamma_moist( ptair(:,:), phumi(:,:) ) * rn_zqt
465	ENDIF
466
467	SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point
468	!
469	CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! NCAR-COREv2
470	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
471	CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! COARE v3.0
472	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
473	CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! COARE v3.5
474	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
475	CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! ECMWF
476	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
477	CASE DEFAULT
478	CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
479	END SELECT
480
481
482	IF ( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp
483	!! FL do we need this multiplication by tmask ... ???
484	DO jj = 1, jpj
485	DO ji = 1, jpi
486	zztmp = zU_zu(ji,jj) * tmask(ji,jj,1)
487	wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod
488	pcd_du(ji,jj) = zztmp * Cd_atm(ji,jj)
489	psen(ji,jj) = zztmp * Ch_atm(ji,jj)
490	pevp(ji,jj) = zztmp * Ce_atm(ji,jj)
491	END DO
492	END DO
493
494	ELSE !== BLK formulation ==! turbulent fluxes computation
495
496	! ! Compute true air density :
497	IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrhogravrn_zu from SLP...)
498	zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , pslp(:,:) )
499	ELSE ! At zt:
500	zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) )
501	END IF
502
503	DO jj = 1, jpj
504	DO ji = 1, jpi
505	zztmp = zrhoa(ji,jj) * zU_zu(ji,jj)
506	!!gm to be done by someone: check the efficiency of the call of cp_air within do loops
507	psen (ji,jj) = cp_air( q_zu(ji,jj) ) * zztmp * Ch_atm(ji,jj) * ( zst(ji,jj) - t_zu(ji,jj) )
508	pevp (ji,jj) = rn_efacMAX( 0._wp, zztmp Ce_atm(ji,jj) * ( pssq(ji,jj) - q_zu(ji,jj) ) )
509	zztmp = zztmp * Cd_atm(ji,jj)
510	pcd_du(ji,jj) = zztmp
511	taum (ji,jj) = zztmp * wndm (ji,jj)
512	zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
513	zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
514	END DO
515	END DO
516
517	CALL iom_put( "taum_oce", taum ) ! output wind stress module
518
519	! ... utau, vtau at U- and V_points, resp.
520	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
521	! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
522	DO jj = 1, jpjm1
523	DO ji = 1, fs_jpim1
524	utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) &
525	& * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
526	vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) &
527	& * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
528	END DO
529	END DO
530	CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
531
532	IF(ln_ctl) THEN
533	CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ')
534	CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, &
535	& tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask )
536	ENDIF
537	!
538	ENDIF
539	!
540	IF(ln_ctl) THEN
541	CALL prt_ctl( tab2d_1=pevp , clinfo1=' blk_oce_1: pevp : ' )
542	CALL prt_ctl( tab2d_1=psen , clinfo1=' blk_oce_1: psen : ' )
543	CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ' )
544	ENDIF
545	!
546	END SUBROUTINE blk_oce_1
547
548
549	SUBROUTINE blk_oce_2( kt, ptair, pqsr, pqlw, pprec, & ! <<= in
550	& psnow, pst , psen, pevp ) ! <<= in
551	!!---------------------------------------------------------------------
552	!! * ROUTINE blk_oce_2 *
553	!!
554	!! ** Purpose : finalize the momentum, heat and freshwater fluxes computation
555	!! at the ocean surface at each time step knowing Cd, Ch, Ce and
556	!! atmospheric variables (from ABL or external data)
557	!!
558	!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
559	!! - vtau : j-component of the stress at V-point (N/m2)
560	!! - taum : Wind stress module at T-point (N/m2)
561	!! - wndm : Wind speed module at T-point (m/s)
562	!! - qsr : Solar heat flux over the ocean (W/m2)
563	!! - qns : Non Solar heat flux over the ocean (W/m2)
564	!! - emp : evaporation minus precipitation (kg/m2/s)
565	!!---------------------------------------------------------------------
566	INTEGER , INTENT(in) :: kt ! time step index
567	REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair
568	REAL(wp), INTENT(in), DIMENSION(:,:) :: pqsr
569	REAL(wp), INTENT(in), DIMENSION(:,:) :: pqlw
570	REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec
571	REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow
572	REAL(wp), INTENT(in), DIMENSION(:,:) :: pst ! surface temperature [Celcius]
573	REAL(wp), INTENT(in), DIMENSION(:,:) :: psen
574	REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp
575	!
576	INTEGER :: ji, jj ! dummy loop indices
577	REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable
578	REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! long wave and sensible heat fluxes
579	REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat fluxes and evaporation
580	REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin
581	!!---------------------------------------------------------------------
582	!
583	! local scalars ( place there for vector optimisation purposes)
584	zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
585
586	! ----------------------------------------------------------------------------- !
587	! I Radiative FLUXES !
588	! ----------------------------------------------------------------------------- !
589
590	! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
591	zztmp = 1._wp - albo
592	IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( pqsr(:,:) ) * tmask(:,:,1)
593	ELSE ; qsr(:,:) = zztmp * pqsr(:,:) * tmask(:,:,1)
594	ENDIF
595
596	zqlw(:,:) = ( pqlw(:,:) - Stef * zst(:,:)zst(:,:)zst(:,:)zst(:,:) ) tmask(:,:,1) ! Long Wave
597
598	! Turbulent fluxes over ocean
599	! -----------------------------
600
601	zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:) ! Latent Heat flux
602
603	IF(ln_ctl) THEN
604	CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ' )
605	CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' )
606	ENDIF
607
608	! ----------------------------------------------------------------------------- !
609	! III Total FLUXES !
610	! ----------------------------------------------------------------------------- !
611	!
612	emp (:,:) = ( pevp(:,:) - pprec(:,:) * rn_pfac ) * tmask(:,:,1) ! mass flux (evap. - precip.)
613	!
614	zz1 = rn_pfac * rLfus
615	zz2 = rn_pfac * rcp
616	zz3 = rn_pfac * rcpi
617	qns(:,:) = ( zqlw(:,:) - psen(:,:) - zqla(:,:) & ! Downward Non Solar
618	& - psnow(:,:) * zztmp & ! remove latent melting heat for solid precip
619	& - pevp(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST
620	& + ( pprec(:,:) - psnow(:,:) ) * ( ptair(:,:) - rt0 ) * zz2 & ! add liquid precip heat content at Tair
621	& + psnow(:,:) * ( MIN( ptair(:,:), rt0 ) - rt0 ) * zz3 & ! add solid precip heat content at min(Tair,Tsnow)
622	& ) * tmask(:,:,1)
623	!
624	#if defined key_si3
625	qns_oce(:,:) = zqlw(:,:) - psen(:,:) - zqla(:,:) ! non solar without emp (only needed by SI3)
626	qsr_oce(:,:) = qsr(:,:)
627	#endif
628	!
629	IF ( nn_ice == 0 ) THEN
630	CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean
631	CALL iom_put( "qsb_oce" , - psen ) ! output downward sensible heat over the ocean
632	CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean
633	CALL iom_put( "qemp_oce", qns-zqlw+psen+zqla ) ! output downward heat content of E-P over the ocean
634	CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean
635	CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean
636	CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean
637	tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
638	sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
639	CALL iom_put( 'snowpre', sprecip ) ! Snow
640	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
641	ENDIF
642	!
643	IF(ln_ctl) THEN
644	CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ')
645	CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ')
646	CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ')
647	ENDIF
648	!
649	END SUBROUTINE blk_oce_2
650
651
652	FUNCTION rho_air( ptak, pqa, pslp )
653	!!-------------------------------------------------------------------------------
654	!! * FUNCTION rho_air *
655	!!
656	!! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
657	!!
658	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
659	!!-------------------------------------------------------------------------------
660	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
661	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
662	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa]
663	REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3]
664	!!-------------------------------------------------------------------------------
665	!
666	rho_air = pslp / ( R_dryptak ( 1._wp + rctv0*pqa ) )
667	!
668	END FUNCTION rho_air
669
670
671	FUNCTION cp_air_0d( pqa )
672	!!-------------------------------------------------------------------------------
673	!! * FUNCTION cp_air *
674	!!
675	!! ** Purpose : Compute specific heat (Cp) of moist air
676	!!
677	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
678	!!-------------------------------------------------------------------------------
679	REAL(wp), INTENT(in) :: pqa ! air specific humidity [kg/kg]
680	REAL(wp) :: cp_air_0d! specific heat of moist air [J/K/kg]
681	!!-------------------------------------------------------------------------------
682	!
683	Cp_air_0d = Cp_dry + Cp_vap * pqa
684	!
685	END FUNCTION cp_air_0d
686
687
688	FUNCTION cp_air_2d( pqa )
689	!!-------------------------------------------------------------------------------
690	!! * FUNCTION cp_air *
691	!!
692	!! ** Purpose : Compute specific heat (Cp) of moist air
693	!!
694	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
695	!!-------------------------------------------------------------------------------
696	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
697	REAL(wp), DIMENSION(jpi,jpj) :: cp_air_2d! specific heat of moist air [J/K/kg]
698	!!-------------------------------------------------------------------------------
699	!
700	Cp_air_2d = Cp_dry + Cp_vap * pqa
701	!
702	END FUNCTION cp_air_2d
703
704
705	FUNCTION q_sat( ptak, pslp )
706	!!----------------------------------------------------------------------------------
707	!! * FUNCTION q_sat *
708	!!
709	!! ** Purpose : Specific humidity at saturation in [kg/kg]
710	!! Based on accurate estimate of "e_sat"
711	!! aka saturation water vapor (Goff, 1957)
712	!!
713	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
714	!!----------------------------------------------------------------------------------
715	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
716	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa]
717	REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg]
718	!
719	INTEGER :: ji, jj ! dummy loop indices
720	REAL(wp) :: ze_sat, ztmp ! local scalar
721	!!----------------------------------------------------------------------------------
722	!
723	DO jj = 1, jpj
724	DO ji = 1, jpi
725	!
726	ztmp = rt0 / ptak(ji,jj)
727	!
728	! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
729	ze_sat = 10.*( 10.79574(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) &
730	& + 1.5047510.(-4)(1. - 10.*(-8.2969(ptak(ji,jj)/rt0 - 1.)) ) &
731	& + 0.4287310.(-3)(10.*(4.76955(1. - ztmp)) - 1.) + 0.78614 )
732	!
733	q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wppslp(ji,jj) - (1._wp - reps0)ze_sat ) ! 0.01 because SLP is in [Pa]
734	!
735	END DO
736	END DO
737	!
738	END FUNCTION q_sat
739
740
741	FUNCTION gamma_moist( ptak, pqa )
742	!!----------------------------------------------------------------------------------
743	!! * FUNCTION gamma_moist *
744	!!
745	!! ** Purpose : Compute the moist adiabatic lapse-rate.
746	!! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
747	!! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
748	!!
749	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
750	!!----------------------------------------------------------------------------------
751	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
752	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg]
753	REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate
754	!
755	INTEGER :: ji, jj ! dummy loop indices
756	REAL(wp) :: zrv, ziRT ! local scalar
757	!!----------------------------------------------------------------------------------
758	!
759	DO jj = 1, jpj
760	DO ji = 1, jpi
761	zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
762	ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT
763	gamma_moist(ji,jj) = grav * ( 1. + rLevapzrvziRT ) / ( Cp_dry + rLevaprLevapzrvreps0ziRT/ptak(ji,jj) )
764	END DO
765	END DO
766	!
767	END FUNCTION gamma_moist
768
769
770	FUNCTION L_vap( psst )
771	!!---------------------------------------------------------------------------------
772	!! * FUNCTION L_vap *
773	!!
774	!! ** Purpose : Compute the latent heat of vaporization of water from temperature
775	!!
776	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
777	!!----------------------------------------------------------------------------------
778	REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg]
779	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]
780	!!----------------------------------------------------------------------------------
781	!
782	L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6
783	!
784	END FUNCTION L_vap
785
786	#if defined key_si3
787	!!----------------------------------------------------------------------
788	!! 'key_si3' SI3 sea-ice model
789	!!----------------------------------------------------------------------
790	!! blk_ice_tau : provide the air-ice stress
791	!! blk_ice_flx : provide the heat and mass fluxes at air-ice interface
792	!! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
793	!! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
794	!! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
795	!!----------------------------------------------------------------------
796
797	SUBROUTINE blk_ice_tau
798	!!---------------------------------------------------------------------
799	!! * ROUTINE blk_ice_tau *
800	!!
801	!! ** Purpose : provide the surface boundary condition over sea-ice
802	!!
803	!! ** Method : compute momentum using bulk formulation
804	!! formulea, ice variables and read atmospheric fields.
805	!! NB: ice drag coefficient is assumed to be a constant
806	!!---------------------------------------------------------------------
807	INTEGER :: ji, jj ! dummy loop indices
808	REAL(wp) :: zwndi_f , zwndj_f, zwnorm_f ! relative wind module and components at F-point
809	REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point
810	REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! transfer coefficient for momentum (tau)
811	!!---------------------------------------------------------------------
812	!
813	! set transfer coefficients to default sea-ice values
814	Cd_atm(:,:) = Cd_ice
815	Ch_atm(:,:) = Cd_ice
816	Ce_atm(:,:) = Cd_ice
817
818	! ------------------------------------------------------------ !
819	! Wind module relative to the moving ice ( U10m - U_ice ) !
820	! ------------------------------------------------------------ !
821	! C-grid ice dynamics : U & V-points (same as ocean)
822	DO jj = 2, jpjm1
823	DO ji = fs_2, fs_jpim1 ! vect. opt.
824	zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) )
825	zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )
826	wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
827	END DO
828	END DO
829	CALL lbc_lnk( 'sbcblk', wndm_ice, 'T', 1. )
830	!
831	! Make ice-atm. drag dependent on ice concentration
832	IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations
833	CALL Cdn10_Lupkes2012( Cd_atm )
834	Ch_atm(:,:) = Cd_atm(:,:) ! momentum and heat transfer coef. are considered identical
835	ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations
836	CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )
837	ENDIF
838
839	!! CALL iom_put( "Cd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef.
840	!! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef.
841
842	! local scalars ( place there for vector optimisation purposes)
843	! Computing density of air! Way denser that 1.2 over sea-ice !!!
844	zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
845
846	! ------------------------------------------------------------ !
847	! Wind stress relative to the moving ice ( U10m - U_ice ) !
848	! ------------------------------------------------------------ !
849	! C-grid ice dynamics : U & V-points (same as ocean)
850	DO jj = 2, jpjm1
851	DO ji = fs_2, fs_jpim1 ! vect. opt.
852	utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) &
853	& * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
854	vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) &
855	& * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
856	END DO
857	END DO
858	CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
859	!
860	!
861	IF(ln_ctl) THEN
862	CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ')
863	CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
864	ENDIF
865	!
866	END SUBROUTINE blk_ice_tau
867
868
869	SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
870	!!---------------------------------------------------------------------
871	!! * ROUTINE blk_ice_flx *
872	!!
873	!! ** Purpose : provide the heat and mass fluxes at air-ice interface
874	!!
875	!! ** Method : compute heat and freshwater exchanged
876	!! between atmosphere and sea-ice using bulk formulation
877	!! formulea, ice variables and read atmmospheric fields.
878	!!
879	!! caution : the net upward water flux has with mm/day unit
880	!!---------------------------------------------------------------------
881	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature
882	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness
883	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness
884	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies)
885	!!
886	INTEGER :: ji, jj, jl ! dummy loop indices
887	REAL(wp) :: zst3 ! local variable
888	REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
889	REAL(wp) :: zztmp, z1_rLsub ! - -
890	REAL(wp) :: zfr1, zfr2 ! local variables
891	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature
892	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice
893	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice
894	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice
895	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice
896	REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3)
897	REAL(wp), DIMENSION(jpi,jpj) :: zrhoa
898	!!---------------------------------------------------------------------
899	!
900	zcoef_dqlw = 4.0 * 0.95 * Stef ! local scalars
901	zcoef_dqla = -Ls * 11637800. * (-5897.8)
902	!
903	zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
904	!
905	zztmp = 1. / ( 1. - albo )
906	WHERE( ptsu(:,:,:) /= 0._wp ) ; z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
907	ELSEWHERE ; z1_st(:,:,:) = 0._wp
908	END WHERE
909	! ! ========================== !
910	DO jl = 1, jpl ! Loop over ice categories !
911	! ! ========================== !
912	DO jj = 1 , jpj
913	DO ji = 1, jpi
914	! ----------------------------!
915	! I Radiative FLUXES !
916	! ----------------------------!
917	zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
918	! Short Wave (sw)
919	qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
920	! Long Wave (lw)
921	z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
922	! lw sensitivity
923	z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
924
925	! ----------------------------!
926	! II Turbulent FLUXES !
927	! ----------------------------!
928
929	! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
930	! Sensible Heat
931	z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))
932	! Latent Heat
933	qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ch_atm(ji,jj) * wndm_ice(ji,jj) * &
934	& ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
935	! Latent heat sensitivity for ice (Dqla/Dt)
936	IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
937	dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) * &
938	& z1_st(ji,jj,jl)z1_st(ji,jj,jl) EXP(-5897.8 * z1_st(ji,jj,jl))
939	ELSE
940	dqla_ice(ji,jj,jl) = 0._wp
941	ENDIF
942
943	! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
944	z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)
945
946	! ----------------------------!
947	! III Total FLUXES !
948	! ----------------------------!
949	! Downward Non Solar flux
950	qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
951	! Total non solar heat flux sensitivity for ice
952	dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
953	END DO
954	!
955	END DO
956	!
957	END DO
958	!
959	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s]
960	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s]
961	CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation
962	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
963
964	! --- evaporation --- !
965	z1_rLsub = 1._wp / rLsub
966	evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation
967	devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT
968	zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean
969
970	! --- evaporation minus precipitation --- !
971	zsnw(:,:) = 0._wp
972	CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing
973	emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
974	emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
975	emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
976
977	! --- heat flux associated with emp --- !
978	qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst
979	& + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & ! liquid precip at Tair
980	& + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow)
981	& ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
982	qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only)
983	& ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
984
985	! --- total solar and non solar fluxes --- !
986	qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) &
987	& + qemp_ice(:,:) + qemp_oce(:,:)
988	qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
989
990	! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
991	qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
992
993	! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
994	DO jl = 1, jpl
995	qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)( ( Tice - rt0 ) rcpi * tmask(:,:,1) )
996	! ! But we do not have Tice => consider it at 0degC => evap=0
997	END DO
998
999	! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
1000	zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) ! transmission when hi>10cm
1001	zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1
1002	!
1003	WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
1004	qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
1005	ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm
1006	qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
1007	ELSEWHERE ! zero when hs>0
1008	qtr_ice_top(:,:,:) = 0._wp
1009	END WHERE
1010	!
1011	IF(ln_ctl) THEN
1012	CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl)
1013	CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
1014	CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl)
1015	CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl)
1016	CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl)
1017	CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ')
1018	ENDIF
1019	!
1020	END SUBROUTINE blk_ice_flx
1021
1022
1023	SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
1024	!!---------------------------------------------------------------------
1025	!! * ROUTINE blk_ice_qcn *
1026	!!
1027	!! ** Purpose : Compute surface temperature and snow/ice conduction flux
1028	!! to force sea ice / snow thermodynamics
1029	!! in the case conduction flux is emulated
1030	!!
1031	!! ** Method : compute surface energy balance assuming neglecting heat storage
1032	!! following the 0-layer Semtner (1976) approach
1033	!!
1034	!! ** Outputs : - ptsu : sea-ice / snow surface temperature (K)
1035	!! - qcn_ice : surface inner conduction flux (W/m2)
1036	!!
1037	!!---------------------------------------------------------------------
1038	LOGICAL , INTENT(in ) :: ld_virtual_itd ! single-category option
1039	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature
1040	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: ptb ! sea ice base temperature
1041	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phs ! snow thickness
1042	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phi ! sea ice thickness
1043	!
1044	INTEGER , PARAMETER :: nit = 10 ! number of iterations
1045	REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction
1046	!
1047	INTEGER :: ji, jj, jl ! dummy loop indices
1048	INTEGER :: iter ! local integer
1049	REAL(wp) :: zfac, zfac2, zfac3 ! local scalars
1050	REAL(wp) :: zkeff_h, ztsu, ztsu0 !
1051	REAL(wp) :: zqc, zqnet !
1052	REAL(wp) :: zhe, zqa0 !
1053	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor
1054	!!---------------------------------------------------------------------
1055
1056	! -------------------------------------!
1057	! I Enhanced conduction factor !
1058	! -------------------------------------!
1059	! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
1060	! Fichefet and Morales Maqueda, JGR 1997
1061	!
1062	zgfac(:,:,:) = 1._wp
1063
1064	IF( ld_virtual_itd ) THEN
1065	!
1066	zfac = 1._wp / ( rn_cnd_s + rcnd_i )
1067	zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
1068	zfac3 = 2._wp / zepsilon
1069	!
1070	DO jl = 1, jpl
1071	DO jj = 1 , jpj
1072	DO ji = 1, jpi
1073	zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness
1074	IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
1075	END DO
1076	END DO
1077	END DO
1078	!
1079	ENDIF
1080
1081	! -------------------------------------------------------------!
1082	! II Surface temperature and conduction flux !
1083	! -------------------------------------------------------------!
1084	!
1085	zfac = rcnd_i * rn_cnd_s
1086	!
1087	DO jl = 1, jpl
1088	DO jj = 1 , jpj
1089	DO ji = 1, jpi
1090	!
1091	zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness
1092	& ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
1093	ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature
1094	ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature
1095	zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
1096	!
1097	DO iter = 1, nit ! --- Iterative loop
1098	zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards)
1099	zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget
1100	ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update
1101	END DO
1102	!
1103	ptsu (ji,jj,jl) = MIN( rt0, ztsu )
1104	qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
1105	qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
1106	qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) ) &
1107	& * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
1108
1109	! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
1110	hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl)
1111
1112	END DO
1113	END DO
1114	!
1115	END DO
1116	!
1117	END SUBROUTINE blk_ice_qcn
1118
1119
1120	SUBROUTINE Cdn10_Lupkes2012( Cd )
1121	!!----------------------------------------------------------------------
1122	!! * ROUTINE Cdn10_Lupkes2012 *
1123	!!
1124	!! ** Purpose : Recompute the neutral air-ice drag referenced at 10m
1125	!! to make it dependent on edges at leads, melt ponds and flows.
1126	!! After some approximations, this can be resumed to a dependency
1127	!! on ice concentration.
1128	!!
1129	!! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50)
1130	!! with the highest level of approximation: level4, eq.(59)
1131	!! The generic drag over a cell partly covered by ice can be re-written as follows:
1132	!!
1133	!! Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)*(nu+1/(10beta)) * A**mu
1134	!!
1135	!! Ce = 2.23e-3 , as suggested by Lupkes (eq. 59)
1136	!! nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
1137	!! A is the concentration of ice minus melt ponds (if any)
1138	!!
1139	!! This new drag has a parabolic shape (as a function of A) starting at
1140	!! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
1141	!! and going down to Cdi(say 1.4e-3) for A=1
1142	!!
1143	!! It is theoretically applicable to all ice conditions (not only MIZ)
1144	!! => see Lupkes et al (2013)
1145	!!
1146	!! ** References : Lupkes et al. JGR 2012 (theory)
1147	!! Lupkes et al. GRL 2013 (application to GCM)
1148	!!
1149	!!----------------------------------------------------------------------
1150	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd
1151	REAL(wp), PARAMETER :: zCe = 2.23e-03_wp
1152	REAL(wp), PARAMETER :: znu = 1._wp
1153	REAL(wp), PARAMETER :: zmu = 1._wp
1154	REAL(wp), PARAMETER :: zbeta = 1._wp
1155	REAL(wp) :: zcoef
1156	!!----------------------------------------------------------------------
1157	zcoef = znu + 1._wp / ( 10._wp * zbeta )
1158
1159	! generic drag over a cell partly covered by ice
1160	!!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag
1161	!! & Cd_ice * at_i_b(:,:) + & ! pure ice drag
1162	!! & zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**zmu ! change due to sea-ice morphology
1163
1164	! ice-atm drag
1165	Cd(:,:) = Cd_ice + & ! pure ice drag
1166	& zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology
1167
1168	END SUBROUTINE Cdn10_Lupkes2012
1169
1170
1171	SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
1172	!!----------------------------------------------------------------------
1173	!! * ROUTINE Cdn10_Lupkes2015 *
1174	!!
1175	!! ** pUrpose : Alternative turbulent transfert coefficients formulation
1176	!! between sea-ice and atmosphere with distinct momentum
1177	!! and heat coefficients depending on sea-ice concentration
1178	!! and atmospheric stability (no meltponds effect for now).
1179	!!
1180	!! ** Method : The parameterization is adapted from Lupkes et al. (2015)
1181	!! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
1182	!! it considers specific skin and form drags (Andreas et al. 2010)
1183	!! to compute neutral transfert coefficients for both heat and
1184	!! momemtum fluxes. Atmospheric stability effect on transfert
1185	!! coefficient is also taken into account following Louis (1979).
1186	!!
1187	!! ** References : Lupkes et al. JGR 2015 (theory)
1188	!! Lupkes et al. ECHAM6 documentation 2015 (implementation)
1189	!!
1190	!!----------------------------------------------------------------------
1191	!
1192	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd
1193	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Ch
1194	REAL(wp), DIMENSION(jpi,jpj) :: ztm_su, zst, zqo_sat, zqi_sat
1195	!
1196	! ECHAM6 constants
1197	REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m]
1198	REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m]
1199	REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m]
1200	REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 41
1201	REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 41
1202	REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 13
1203	REAL(wp), PARAMETER :: zc2 = zc * zc
1204	REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 14
1205	REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 30
1206	REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 51
1207	REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 56
1208	REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 26
1209	REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 26
1210	REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma
1211	REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp
1212	!
1213	INTEGER :: ji, jj ! dummy loop indices
1214	REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu
1215	REAL(wp) :: zrib_o, zrib_i
1216	REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice
1217	REAL(wp) :: zChn_skin_ice, zChn_form_ice
1218	REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw
1219	REAL(wp) :: zCdn_form_tmp
1220	!!----------------------------------------------------------------------
1221
1222	! mean temperature
1223	WHERE( at_i_b(:,:) > 1.e-20 ) ; ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
1224	ELSEWHERE ; ztm_su(:,:) = rt0
1225	ENDWHERE
1226
1227	! Momentum Neutral Transfert Coefficients (should be a constant)
1228	zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40
1229	zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7
1230	zCdn_ice = zCdn_skin_ice ! Eq. 7 (cf Lupkes email for details)
1231	!zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32)
1232
1233	! Heat Neutral Transfert Coefficients
1234	zChn_skin_ice = vkarmn*2 / ( LOG( rn_zu / z0_ice + 1._wp ) LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) ) ! Eq. 50 + Eq. 52 (cf Lupkes email for details)
1235
1236	! Atmospheric and Surface Variables
1237	zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin
1238	zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ocean [kg/kg]
1239	zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg]
1240	!
1241	DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility
1242	DO ji = fs_2, fs_jpim1
1243	! Virtual potential temperature [K]
1244	zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean
1245	zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice
1246	zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu
1247
1248	! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
1249	zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean
1250	zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice
1251
1252	! Momentum and Heat Neutral Transfert Coefficients
1253	zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40
1254	zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53
1255
1256	! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
1257	z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
1258	z0i = z0_skin_ice ! over ice (cf Lupkes email for details)
1259	IF( zrib_o <= 0._wp ) THEN
1260	zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) ) ! Eq. 10
1261	zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )*r1_3 / ( Chn_oce(ji,jj) MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26
1262	& )zgamma )z1_gamma
1263	ELSE
1264	zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12
1265	zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28
1266	ENDIF
1267
1268	IF( zrib_i <= 0._wp ) THEN
1269	zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9
1270	zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25
1271	ELSE
1272	zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11
1273	zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27
1274	ENDIF
1275
1276	! Momentum Transfert Coefficients (Eq. 38)
1277	Cd(ji,jj) = zCdn_skin_ice * zfmi + &
1278	& zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
1279
1280	! Heat Transfert Coefficients (Eq. 49)
1281	Ch(ji,jj) = zChn_skin_ice * zfhi + &
1282	& zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
1283	!
1284	END DO
1285	END DO
1286	CALL lbc_lnk_multi( 'sbcblk', Cd, 'T', 1., Ch, 'T', 1. )
1287	!
1288	END SUBROUTINE Cdn10_Lupkes2015
1289
1290	#endif
1291
1292	!!======================================================================
1293	END MODULE sbcblk

Note: See TracBrowser for help on using the repository browser.

Download in other formats: