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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sbcblk.F90 in NEMO/trunk/cfgs/SPITZ12/MY_SRC – NEMO

Context Navigation

source: NEMO/trunk/cfgs/SPITZ12/MY_SRC/sbcblk.F90 @ 9598

Last change on this file since 9598 was 9598, checked in by nicolasmartin, 6 years ago
Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing Standardisation of routines header with version 4.0 and year 2018 Fix for some broken symlinks
File size: 73.0 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_si3
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 lbclnk ! ocean lateral boundary conditions (or mpp link)
60	USE prtctl ! Print control
61
62	IMPLICIT NONE
63	PRIVATE
64
65	PUBLIC sbc_blk_init ! called in sbcmod
66	PUBLIC sbc_blk ! called in sbcmod
67	#if defined key_si3
68	PUBLIC blk_ice_tau ! routine called in iceforcing
69	PUBLIC blk_ice_flx ! routine called in iceforcing
70	PUBLIC blk_ice_qcn ! routine called in iceforcing
71	#endif
72
73	!!Lolo: should ultimately be moved in the module with all physical constants ?
74	!!gm : In principle, yes.
75	REAL(wp), PARAMETER :: Cp_dry = 1005.0 !: Specic heat of dry air, constant pressure [J/K/kg]
76	REAL(wp), PARAMETER :: Cp_vap = 1860.0 !: Specic heat of water vapor, constant pressure [J/K/kg]
77	REAL(wp), PARAMETER :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg]
78	REAL(wp), PARAMETER :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg]
79	REAL(wp), PARAMETER :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
80	REAL(wp), PARAMETER :: rctv0 = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
81
82	!!clem INTEGER , PARAMETER :: jpfld =10 ! maximum number of files to read
83	INTEGER , PARAMETER :: jpfld = 8 ! 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 (./LICENSE)
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	!!clem 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	!!clem jfld = jpfld - COUNT( (/.NOT.lhftau/) )
221	jfld = jpfld
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	IF( slf_i(ifpr)%nfreqh > 0. .AND. MOD( 3600. * slf_i(ifpr)%nfreqh , REAL(nn_fsbc) * rdt) /= 0. ) &
230	& CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', &
231	& ' This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
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_tauwoc .OR. ln_stcor ) ) THEN
240	CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=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_tauwoc .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_tauwoc=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	! local scalars ( place there for vector optimisation purposes)
387	zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
388
389	! ----------------------------------------------------------------------------- !
390	! 0 Wind components and module at T-point relative to the moving ocean !
391	! ----------------------------------------------------------------------------- !
392
393	! ... components ( U10m - U_oce ) at T-point (unmasked)
394	!!gm move zwnd_i (_j) set to zero inside the key_cyclone ???
395	zwnd_i(:,:) = 0._wp
396	zwnd_j(:,:) = 0._wp
397	#if defined key_cyclone
398	CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012)
399	DO jj = 2, jpjm1
400	DO ji = fs_2, fs_jpim1 ! vect. opt.
401	sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
402	sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
403	END DO
404	END DO
405	#endif
406	DO jj = 2, jpjm1
407	DO ji = fs_2, fs_jpim1 ! vect. opt.
408	zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) )
409	zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) )
410	END DO
411	END DO
412	CALL lbc_lnk_multi( zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
413	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
414	wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) &
415	& + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1)
416
417	! ----------------------------------------------------------------------------- !
418	! I Radiative FLUXES !
419	! ----------------------------------------------------------------------------- !
420
421	! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
422	zztmp = 1. - albo
423	IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
424	ELSE ; qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1)
425	ENDIF
426
427	zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)zst(:,:)zst(:,:)zst(:,:) ) tmask(:,:,1) ! Long Wave
428
429	! ----------------------------------------------------------------------------- !
430	! II Turbulent FLUXES !
431	! ----------------------------------------------------------------------------- !
432
433	! ... specific humidity at SST and IST tmask(
434	!!clem zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
435	zsq(:,:) = 0.98 * q_sat( zst(:,:), 101300. )
436	!!
437	!! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
438	!! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
439	!! (since reanalysis products provide T at z, not theta !)
440	ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt
441
442	SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point
443	!
444	CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2
445	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
446	CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0
447	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
448	CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5
449	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
450	CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF
451	& Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
452	CASE DEFAULT
453	CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
454	END SELECT
455
456	! ! Compute true air density :
457	IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrhogravrn_zu from SLP...)
458	!!clem zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) )
459	zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , 101300. )
460	ELSE ! At zt:
461	!!clem zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
462	zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), 101300. )
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	!!clem 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_multi( utau, 'U', -1., vtau, 'V', -1. )
494
495	! Turbulent fluxes over ocean
496	! -----------------------------
497
498	! zqla used as temporary array, for rho*U (common term of bulk formulae):
499	zqla(:,:) = zrhoa(:,:) * zU_zu(:,:)
500
501	IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
502	!! q_air and t_air are given at 10m (wind reference height)
503	zevap(:,:) = rn_efacMAX( 0._wp, zqla(:,:)Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed
504	zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))zqla(:,:)Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed
505	ELSE
506	!! q_air and t_air are not given at 10m (wind reference height)
507	! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
508	zevap(:,:) = rn_efacMAX( 0._wp, zqla(:,:)Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed
509	zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))zqla(:,:)Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed
510	ENDIF
511
512	zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:) ! Latent Heat flux
513
514
515	IF(ln_ctl) THEN
516	CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce : ' )
517	CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce : ' )
518	CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' )
519	CALL prt_ctl( tab2d_1=zsq , clinfo1=' blk_oce: zsq : ', tab2d_2=zst, clinfo2=' zst : ' )
520	CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, &
521	& tab2d_2=vtau , clinfo2= ' vtau : ', mask2=vmask )
522	CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce: wndm : ')
523	CALL prt_ctl( tab2d_1=zst , clinfo1=' blk_oce: zst : ')
524	ENDIF
525
526	! ----------------------------------------------------------------------------- !
527	! III Total FLUXES !
528	! ----------------------------------------------------------------------------- !
529	!
530	emp (:,:) = ( zevap(:,:) & ! mass flux (evap. - precip.)
531	& - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1)
532	!
533	qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar
534	& - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus & ! remove latent melting heat for solid precip
535	& - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST
536	& + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac & ! add liquid precip heat content at Tair
537	& * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp &
538	& + sf(jp_snow)%fnow(:,:,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)
539	& * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1)
540	!
541	#if defined key_si3
542	qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) ! non solar without emp (only needed by LIM3)
543	qsr_oce(:,:) = qsr(:,:)
544	#endif
545	!
546	IF ( nn_ice == 0 ) THEN
547	CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean
548	CALL iom_put( "qsb_oce" , - zqsb ) ! output downward sensible heat over the ocean
549	CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean
550	CALL iom_put( "qemp_oce", qns-zqlw+zqsb+zqla ) ! output downward heat content of E-P over the ocean
551	CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean
552	CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean
553	CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean
554	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac ! output total precipitation [kg/m2/s]
555	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac ! output solid precipitation [kg/m2/s]
556	CALL iom_put( 'snowpre', sprecip ) ! Snow
557	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
558	ENDIF
559	!
560	IF(ln_ctl) THEN
561	CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=zqlw , clinfo2=' zqlw : ')
562	CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ')
563	CALL prt_ctl(tab2d_1=pst , clinfo1=' blk_oce: pst : ', tab2d_2=emp , clinfo2=' emp : ')
564	CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, &
565	& tab2d_2=vtau , clinfo2= ' vtau : ' , mask2=vmask )
566	ENDIF
567	!
568	END SUBROUTINE blk_oce
569
570
571
572	FUNCTION rho_air( ptak, pqa, pslp )
573	!!-------------------------------------------------------------------------------
574	!! * FUNCTION rho_air *
575	!!
576	!! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
577	!!
578	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
579	!!-------------------------------------------------------------------------------
580	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
581	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
582	!!clem REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa]
583	REAL(wp), INTENT(in) :: pslp ! pressure in [Pa]
584	REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3]
585	!!-------------------------------------------------------------------------------
586	!
587	rho_air = pslp / ( R_dryptak ( 1._wp + rctv0*pqa ) )
588	!
589	END FUNCTION rho_air
590
591
592	FUNCTION cp_air( pqa )
593	!!-------------------------------------------------------------------------------
594	!! * FUNCTION cp_air *
595	!!
596	!! ** Purpose : Compute specific heat (Cp) of moist air
597	!!
598	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
599	!!-------------------------------------------------------------------------------
600	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
601	REAL(wp), DIMENSION(jpi,jpj) :: cp_air ! specific heat of moist air [J/K/kg]
602	!!-------------------------------------------------------------------------------
603	!
604	Cp_air = Cp_dry + Cp_vap * pqa
605	!
606	END FUNCTION cp_air
607
608
609	FUNCTION q_sat( ptak, pslp )
610	!!----------------------------------------------------------------------------------
611	!! * FUNCTION q_sat *
612	!!
613	!! ** Purpose : Specific humidity at saturation in [kg/kg]
614	!! Based on accurate estimate of "e_sat"
615	!! aka saturation water vapor (Goff, 1957)
616	!!
617	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
618	!!----------------------------------------------------------------------------------
619	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
620	!!clem REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa]
621	REAL(wp), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa]
622	REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg]
623	!
624	INTEGER :: ji, jj ! dummy loop indices
625	REAL(wp) :: ze_sat, ztmp ! local scalar
626	!!----------------------------------------------------------------------------------
627	!
628	DO jj = 1, jpj
629	DO ji = 1, jpi
630	!
631	ztmp = rt0 / ptak(ji,jj)
632	!
633	! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
634	ze_sat = 10.*( 10.79574(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) &
635	& + 1.5047510.(-4)(1. - 10.*(-8.2969(ptak(ji,jj)/rt0 - 1.)) ) &
636	& + 0.4287310.(-3)(10.*(4.76955(1. - ztmp)) - 1.) + 0.78614 )
637	!
638	!!clem q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wppslp(ji,jj) - (1._wp - reps0)ze_sat ) ! 0.01 because SLP is in [Pa]
639	q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wppslp - (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_si3
693	!!----------------------------------------------------------------------
694	!! 'key_si3' 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	! set transfer coefficients to default sea-ice values
720	Cd_atm(:,:) = Cd_ice
721	Ch_atm(:,:) = Cd_ice
722	Ce_atm(:,:) = Cd_ice
723
724	wndm_ice(:,:) = 0._wp !!gm brutal....
725
726	! ------------------------------------------------------------ !
727	! Wind module relative to the moving ice ( U10m - U_ice ) !
728	! ------------------------------------------------------------ !
729	SELECT CASE( cp_ice_msh )
730	CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation)
731	! and scalar wind at T-point ( = \| U10m - U_ice \| ) (masked)
732	DO jj = 2, jpjm1
733	DO ji = 2, jpim1 ! B grid : NO vector opt
734	! ... scalar wind at T-point (fld being at T-point)
735	zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( u_ice(ji,jj+1) + u_ice(ji+1,jj+1) &
736	& + u_ice(ji,jj ) + u_ice(ji+1,jj ) )
737	zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) &
738	& + v_ice(ji,jj ) + v_ice(ji+1,jj ) )
739	wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
740	END DO
741	END DO
742	CALL lbc_lnk( wndm_ice, 'T', 1. )
743	!
744	CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean)
745	DO jj = 2, jpjm1
746	DO ji = fs_2, fs_jpim1 ! vect. opt.
747	zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) )
748	zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )
749	wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
750	END DO
751	END DO
752	CALL lbc_lnk( wndm_ice, 'T', 1. )
753	!
754	END SELECT
755
756	! Make ice-atm. drag dependent on ice concentration
757	IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations
758	CALL Cdn10_Lupkes2012( Cd_atm )
759	Ch_atm(:,:) = Cd_atm(:,:) ! momentum and heat transfer coef. are considered identical
760	ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations
761	CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )
762	ENDIF
763
764	!! CALL iom_put( "Cd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef.
765	!! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef.
766
767	! local scalars ( place there for vector optimisation purposes)
768	! Computing density of air! Way denser that 1.2 over sea-ice !!!
769	!!clem zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
770	zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), 101300.)
771
772	!!gm brutal....
773	utau_ice (:,:) = 0._wp
774	vtau_ice (:,:) = 0._wp
775	!!gm end
776
777	! ------------------------------------------------------------ !
778	! Wind stress relative to the moving ice ( U10m - U_ice ) !
779	! ------------------------------------------------------------ !
780	SELECT CASE( cp_ice_msh )
781	CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation)
782	DO jj = 2, jpjm1
783	DO ji = 2, jpim1 ! B grid : NO vector opt
784	! ... scalar wind at I-point (fld being at T-point)
785	zwndi_f = 0.25 * ( sf(jp_wndi)%fnow(ji-1,jj ,1) + sf(jp_wndi)%fnow(ji ,jj ,1) &
786	& + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji ,jj-1,1) ) - rn_vfac * u_ice(ji,jj)
787	zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) &
788	& + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * v_ice(ji,jj)
789	! ... ice stress at I-point
790	zwnorm_f = zrhoa(ji,jj) * Cd_atm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
791	utau_ice(ji,jj) = zwnorm_f * zwndi_f
792	vtau_ice(ji,jj) = zwnorm_f * zwndj_f
793	END DO
794	END DO
795	CALL lbc_lnk_multi( utau_ice, 'I', -1., vtau_ice, 'I', -1. )
796	!
797	CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean)
798	DO jj = 2, jpjm1
799	DO ji = fs_2, fs_jpim1 ! vect. opt.
800	utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) &
801	& * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
802	vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) &
803	& * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
804	END DO
805	END DO
806	CALL lbc_lnk_multi( utau_ice, 'U', -1., vtau_ice, 'V', -1. )
807	!
808	END SELECT
809	!
810	IF(ln_ctl) THEN
811	CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ')
812	CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
813	ENDIF
814	!
815	END SUBROUTINE blk_ice_tau
816
817
818	SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
819	!!---------------------------------------------------------------------
820	!! * ROUTINE blk_ice_flx *
821	!!
822	!! ** Purpose : provide the heat and mass fluxes at air-ice interface
823	!!
824	!! ** Method : compute heat and freshwater exchanged
825	!! between atmosphere and sea-ice using bulk formulation
826	!! formulea, ice variables and read atmmospheric fields.
827	!!
828	!! caution : the net upward water flux has with mm/day unit
829	!!---------------------------------------------------------------------
830	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature
831	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness
832	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness
833	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies)
834	!!
835	INTEGER :: ji, jj, jl ! dummy loop indices
836	REAL(wp) :: zst3 ! local variable
837	REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
838	REAL(wp) :: zztmp, z1_lsub ! - -
839	REAL(wp) :: zfr1, zfr2 ! local variables
840	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature
841	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice
842	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice
843	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice
844	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice
845	REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3)
846	REAL(wp), DIMENSION(jpi,jpj) :: zrhoa
847	!!---------------------------------------------------------------------
848	!
849	zcoef_dqlw = 4.0 * 0.95 * Stef ! local scalars
850	zcoef_dqla = -Ls * 11637800. * (-5897.8)
851	!
852	!!clem zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
853	zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), 101300. )
854	!
855	zztmp = 1. / ( 1. - albo )
856	WHERE( ptsu(:,:,:) /= 0._wp ) ; z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
857	ELSEWHERE ; z1_st(:,:,:) = 0._wp
858	END WHERE
859	! ! ========================== !
860	DO jl = 1, jpl ! Loop over ice categories !
861	! ! ========================== !
862	DO jj = 1 , jpj
863	DO ji = 1, jpi
864	! ----------------------------!
865	! I Radiative FLUXES !
866	! ----------------------------!
867	zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
868	! Short Wave (sw)
869	qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
870	! Long Wave (lw)
871	z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
872	! lw sensitivity
873	z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
874
875	! ----------------------------!
876	! II Turbulent FLUXES !
877	! ----------------------------!
878
879	! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
880	! Sensible Heat
881	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))
882	! Latent Heat
883	qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ch_atm(ji,jj) * wndm_ice(ji,jj) * &
884	& ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
885	! Latent heat sensitivity for ice (Dqla/Dt)
886	IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
887	dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) * &
888	& z1_st(ji,jj,jl)z1_st(ji,jj,jl) EXP(-5897.8 * z1_st(ji,jj,jl))
889	ELSE
890	dqla_ice(ji,jj,jl) = 0._wp
891	ENDIF
892
893	! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
894	z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)
895
896	! ----------------------------!
897	! III Total FLUXES !
898	! ----------------------------!
899	! Downward Non Solar flux
900	qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
901	! Total non solar heat flux sensitivity for ice
902	dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
903	END DO
904	!
905	END DO
906	!
907	END DO
908	!
909	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac ! total precipitation [kg/m2/s]
910	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac ! solid precipitation [kg/m2/s]
911	CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation
912	CALL iom_put( 'precip' , tprecip ) ! Total precipitation
913
914	! --- evaporation --- !
915	z1_lsub = 1._wp / Lsub
916	evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub ! sublimation
917	devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub ! d(sublimation)/dT
918	zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean
919
920	! --- evaporation minus precipitation --- !
921	zsnw(:,:) = 0._wp
922	CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing
923	emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
924	emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
925	emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
926
927	! --- heat flux associated with emp --- !
928	qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst
929	& + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & ! liquid precip at Tair
930	& + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow)
931	& ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
932	qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only)
933	& ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
934
935	! --- total solar and non solar fluxes --- !
936	qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) &
937	& + qemp_ice(:,:) + qemp_oce(:,:)
938	qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
939
940	! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
941	qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
942
943	! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
944	DO jl = 1, jpl
945	qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)( ( Tice - rt0 ) cpic * tmask(:,:,1) )
946	! ! But we do not have Tice => consider it at 0degC => evap=0
947	END DO
948
949	! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
950	zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) ! transmission when hi>10cm
951	zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1
952	!
953	WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
954	qsr_ice_tr(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
955	ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm
956	qsr_ice_tr(:,:,:) = qsr_ice(:,:,:) * zfr1
957	ELSEWHERE ! zero when hs>0
958	qsr_ice_tr(:,:,:) = 0._wp
959	END WHERE
960	!
961	IF(ln_ctl) THEN
962	CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl)
963	CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
964	CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl)
965	CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl)
966	CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl)
967	CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ')
968	ENDIF
969	!
970	END SUBROUTINE blk_ice_flx
971
972
973	SUBROUTINE blk_ice_qcn( k_virtual_itd, ptsu, ptb, phs, phi )
974	!!---------------------------------------------------------------------
975	!! * ROUTINE blk_ice_qcn *
976	!!
977	!! ** Purpose : Compute surface temperature and snow/ice conduction flux
978	!! to force sea ice / snow thermodynamics
979	!! in the case JULES coupler is emulated
980	!!
981	!! ** Method : compute surface energy balance assuming neglecting heat storage
982	!! following the 0-layer Semtner (1976) approach
983	!!
984	!! ** Outputs : - ptsu : sea-ice / snow surface temperature (K)
985	!! - qcn_ice : surface inner conduction flux (W/m2)
986	!!
987	!!---------------------------------------------------------------------
988	INTEGER , INTENT(in ) :: k_virtual_itd ! single-category option
989	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature
990	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: ptb ! sea ice base temperature
991	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phs ! snow thickness
992	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phi ! sea ice thickness
993	!
994	INTEGER , PARAMETER :: nit = 10 ! number of iterations
995	REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction
996	!
997	INTEGER :: ji, jj, jl ! dummy loop indices
998	INTEGER :: iter ! local integer
999	REAL(wp) :: zfac, zfac2, zfac3 ! local scalars
1000	REAL(wp) :: zkeff_h, ztsu, ztsu0 !
1001	REAL(wp) :: zqc, zqnet !
1002	REAL(wp) :: zhe, zqa0 !
1003	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor
1004	!!---------------------------------------------------------------------
1005
1006	! -------------------------------------!
1007	! I Enhanced conduction factor !
1008	! -------------------------------------!
1009	! Emulates the enhancement of conduction by unresolved thin ice (k_virtual_itd = 1/2)
1010	! Fichefet and Morales Maqueda, JGR 1997
1011	!
1012	zgfac(:,:,:) = 1._wp
1013
1014	SELECT CASE ( k_virtual_itd )
1015	!
1016	CASE ( 1 , 2 )
1017	!
1018	zfac = 1._wp / ( rn_cnd_s + rcdic )
1019	zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
1020	zfac3 = 2._wp / zepsilon
1021	!
1022	DO jl = 1, jpl
1023	DO jj = 1 , jpj
1024	DO ji = 1, jpi
1025	zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac ! Effective thickness
1026	IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
1027	END DO
1028	END DO
1029	END DO
1030	!
1031	END SELECT
1032
1033	! -------------------------------------------------------------!
1034	! II Surface temperature and conduction flux !
1035	! -------------------------------------------------------------!
1036	!
1037	zfac = rcdic * rn_cnd_s
1038	!
1039	DO jl = 1, jpl
1040	DO jj = 1 , jpj
1041	DO ji = 1, jpi
1042	!
1043	zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness
1044	& ( rcdic * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
1045	ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature
1046	ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature
1047	zqa0 = qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
1048	!
1049	DO iter = 1, nit ! --- Iterative loop
1050	zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards)
1051	zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget
1052	ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update
1053	END DO
1054	!
1055	ptsu (ji,jj,jl) = MIN( rt0, ztsu )
1056	qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
1057	qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
1058	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) ) &
1059	& * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
1060
1061	END DO
1062	END DO
1063	!
1064	END DO
1065	!
1066	END SUBROUTINE blk_ice_qcn
1067
1068
1069	SUBROUTINE Cdn10_Lupkes2012( Cd )
1070	!!----------------------------------------------------------------------
1071	!! * ROUTINE Cdn10_Lupkes2012 *
1072	!!
1073	!! ** Purpose : Recompute the neutral air-ice drag referenced at 10m
1074	!! to make it dependent on edges at leads, melt ponds and flows.
1075	!! After some approximations, this can be resumed to a dependency
1076	!! on ice concentration.
1077	!!
1078	!! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50)
1079	!! with the highest level of approximation: level4, eq.(59)
1080	!! The generic drag over a cell partly covered by ice can be re-written as follows:
1081	!!
1082	!! Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)*(nu+1/(10beta)) * A**mu
1083	!!
1084	!! Ce = 2.23e-3 , as suggested by Lupkes (eq. 59)
1085	!! nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
1086	!! A is the concentration of ice minus melt ponds (if any)
1087	!!
1088	!! This new drag has a parabolic shape (as a function of A) starting at
1089	!! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
1090	!! and going down to Cdi(say 1.4e-3) for A=1
1091	!!
1092	!! It is theoretically applicable to all ice conditions (not only MIZ)
1093	!! => see Lupkes et al (2013)
1094	!!
1095	!! ** References : Lupkes et al. JGR 2012 (theory)
1096	!! Lupkes et al. GRL 2013 (application to GCM)
1097	!!
1098	!!----------------------------------------------------------------------
1099	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd
1100	REAL(wp), PARAMETER :: zCe = 2.23e-03_wp
1101	REAL(wp), PARAMETER :: znu = 1._wp
1102	REAL(wp), PARAMETER :: zmu = 1._wp
1103	REAL(wp), PARAMETER :: zbeta = 1._wp
1104	REAL(wp) :: zcoef
1105	!!----------------------------------------------------------------------
1106	zcoef = znu + 1._wp / ( 10._wp * zbeta )
1107
1108	! generic drag over a cell partly covered by ice
1109	!!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag
1110	!! & Cd_ice * at_i_b(:,:) + & ! pure ice drag
1111	!! & zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**zmu ! change due to sea-ice morphology
1112
1113	! ice-atm drag
1114	Cd(:,:) = Cd_ice + & ! pure ice drag
1115	& zCe * ( 1._wp - at_i_b(:,:) )*zcoef at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology
1116
1117	END SUBROUTINE Cdn10_Lupkes2012
1118
1119
1120	SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
1121	!!----------------------------------------------------------------------
1122	!! * ROUTINE Cdn10_Lupkes2015 *
1123	!!
1124	!! ** pUrpose : Alternative turbulent transfert coefficients formulation
1125	!! between sea-ice and atmosphere with distinct momentum
1126	!! and heat coefficients depending on sea-ice concentration
1127	!! and atmospheric stability (no meltponds effect for now).
1128	!!
1129	!! ** Method : The parameterization is adapted from Lupkes et al. (2015)
1130	!! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
1131	!! it considers specific skin and form drags (Andreas et al. 2010)
1132	!! to compute neutral transfert coefficients for both heat and
1133	!! momemtum fluxes. Atmospheric stability effect on transfert
1134	!! coefficient is also taken into account following Louis (1979).
1135	!!
1136	!! ** References : Lupkes et al. JGR 2015 (theory)
1137	!! Lupkes et al. ECHAM6 documentation 2015 (implementation)
1138	!!
1139	!!----------------------------------------------------------------------
1140	!
1141	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd
1142	REAL(wp), DIMENSION(:,:), INTENT(inout) :: Ch
1143	REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat
1144	!
1145	! ECHAM6 constants
1146	REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m]
1147	REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m]
1148	REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m]
1149	REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 41
1150	REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 41
1151	REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 13
1152	REAL(wp), PARAMETER :: zc2 = zc * zc
1153	REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 14
1154	REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 30
1155	REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 51
1156	REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 56
1157	REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 26
1158	REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 26
1159	REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma
1160	REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp
1161	!
1162	INTEGER :: ji, jj ! dummy loop indices
1163	REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu
1164	REAL(wp) :: zrib_o, zrib_i
1165	REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice
1166	REAL(wp) :: zChn_skin_ice, zChn_form_ice
1167	REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw
1168	REAL(wp) :: zCdn_form_tmp
1169	!!----------------------------------------------------------------------
1170
1171	! Momentum Neutral Transfert Coefficients (should be a constant)
1172	zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40
1173	zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7
1174	zCdn_ice = zCdn_skin_ice ! Eq. 7 (cf Lupkes email for details)
1175	!zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32)
1176
1177	! Heat Neutral Transfert Coefficients
1178	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)
1179
1180	! Atmospheric and Surface Variables
1181	zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin
1182	!!clem zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ocean [kg/kg]
1183	!!clem zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg]
1184	zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , 101300. ) ! saturation humidity over ocean [kg/kg]
1185	zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), 101300. ) ! saturation humidity over ice [kg/kg]
1186	!
1187	DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility
1188	DO ji = fs_2, fs_jpim1
1189	! Virtual potential temperature [K]
1190	zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean
1191	zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice
1192	zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu
1193
1194	! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
1195	zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean
1196	zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice
1197
1198	! Momentum and Heat Neutral Transfert Coefficients
1199	zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40
1200	zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53
1201
1202	! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
1203	z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
1204	z0i = z0_skin_ice ! over ice (cf Lupkes email for details)
1205	IF( zrib_o <= 0._wp ) THEN
1206	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
1207	zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )*r1_3 / ( Chn_oce(ji,jj) MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26
1208	& )zgamma )z1_gamma
1209	ELSE
1210	zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12
1211	zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28
1212	ENDIF
1213
1214	IF( zrib_i <= 0._wp ) THEN
1215	zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9
1216	zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25
1217	ELSE
1218	zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11
1219	zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27
1220	ENDIF
1221
1222	! Momentum Transfert Coefficients (Eq. 38)
1223	Cd(ji,jj) = zCdn_skin_ice * zfmi + &
1224	& 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) )
1225
1226	! Heat Transfert Coefficients (Eq. 49)
1227	Ch(ji,jj) = zChn_skin_ice * zfhi + &
1228	& 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) )
1229	!
1230	END DO
1231	END DO
1232	CALL lbc_lnk_multi( Cd, 'T', 1., Ch, 'T', 1. )
1233	!
1234	END SUBROUTINE Cdn10_Lupkes2015
1235
1236	#endif
1237
1238	!!======================================================================
1239	END MODULE sbcblk

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

Download in other formats: