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

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

sbcblk.F90 in NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90 @ 12043

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

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

Download in other formats: