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_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90 @ 11996

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

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

Download in other formats: