Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sbcblk.F90 @ 12324

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/SBC/sbcblk.F90 @ 12324

Download in other formats: