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

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

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

Download in other formats: