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

Last change on this file since 12162 was 12162, checked in by laurent, 4 years ago
Catching up with "dev_r12072_MERGE_OPTION2_2019"
Property svn:keywords set to `Id`
File size: 70.8 KB

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

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

Download in other formats: