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 branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90 @ 8962

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

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

Download in other formats: