Context Navigation

sbcblk.F90 @ 9015

Last change on this file since 9015 was 9015, checked in by acc, 7 years ago
Branch dev_CNRS_2017. Tidy up timing calls. Correct a few inconsistent labels and fix incorrect call in trc_sub_stp. Remove timing from icethd_pnd which is not always called by all processors
File size: 72.1 KB

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

Download in other formats: