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_ASINTER-01-05_merged/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk.F90 @ 12021

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

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

Download in other formats: