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

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

sbcblk.F90 in NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk.F90 @ 12785

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

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

Download in other formats: