Context Navigation

sbcblk.F90 @ 8813

Last change on this file since 8813 was 8813, checked in by gm, 6 years ago
#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8787
File size: 72.5 KB

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

Download in other formats: