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

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

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

Download in other formats: