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_r11643_ENHANCE-11_CEthe_Shaconemo_diags/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11643_ENHANCE-11_CEthe_Shaconemo_diags/src/OCE/SBC/sbcblk.F90 @ 12124

Last change on this file since 12124 was 12124, checked in by clem, 4 years ago
cosmetics
Property svn:keywords set to `Id`
File size: 71.6 KB

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

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

Download in other formats: