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

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

sbcblk.F90 in branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

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

Last change on this file since 8882 was 8882, checked in by flavoni, 6 years ago
dev_CNRS_2017 branch: merged dev_r7881_ENHANCE09_RK3 with trunk r8864
File size: 72.8 KB

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

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

Download in other formats: