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/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/SBC/sbcblk.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 5 years ago
Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
Property svn:keywords set to `Id`
File size: 70.0 KB

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

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

Download in other formats: