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 @ 14717

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

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

Download in other formats: