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/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbcblk.F90 @ 12991

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

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

Download in other formats: