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/releases/r4.0/r4.0-HEAD/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcblk.F90 @ 15613

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

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

Download in other formats: