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_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90 @ 11334

Last change on this file since 11334 was 11334, checked in by flemarie, 5 years ago

Follow-on to the ABL implementation (see ticket #2131)

Modification of diawri.F90 to output ABL variables with IOIPSL (i.e. without key_iomput)
Modification of ice_var_agg in icevar.F90 to allow the use of tm_su in ABL and bulk
Error handling in case ln_abl = .true. and ABL sources were not compiled
ABL now works with SI3 considering an average over ice categories
Update reference namelist to avoid runtime warnings due to nam_tide

Sanity checks with ORCA2 for the following configurations :

1 - ABL src + ln_blk = .true.
2 - ABL src + ln_abl = .true.
3 - no ABL src + ln_blk = .true.

All configurations are MPP reproducible and configurations 1 and 3 results are identical

Property svn:keywords set to Id

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

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

Download in other formats: