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

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

sbcblk.F90 in branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90 @ 7163

Last change on this file since 7163 was 7163, checked in by gm, 7 years ago
#1751 - branch SIMPLIF_6_aerobulk: update option control in sbcmod + uniformization of print in ocean_output (many module involved)
File size: 51.7 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	!!----------------------------------------------------------------------
20
21	!!----------------------------------------------------------------------
22	!! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition
23	!! sbc_blk : bulk formulation as ocean surface boundary condition
24	!! blk_oce : computes momentum, heat and freshwater fluxes over ocean
25	!! blk_ice : computes momentum, heat and freshwater fluxes over sea ice
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	!!----------------------------------------------------------------------
31	USE oce ! ocean dynamics and tracers
32	USE dom_oce ! ocean space and time domain
33	USE phycst ! physical constants
34	USE fldread ! read input fields
35	USE sbc_oce ! Surface boundary condition: ocean fields
36	USE cyclone ! Cyclone 10m wind form trac of cyclone centres
37	USE sbcdcy ! surface boundary condition: diurnal cycle
38	USE sbcwave , ONLY : cdn_wave ! wave module
39	USE sbc_ice ! Surface boundary condition: ice fields
40	USE lib_fortran ! to use key_nosignedzero
41	#if defined key_lim3
42	USE ice , ONLY : u_ice, v_ice, jpl, pfrld, a_i_b
43	USE limthd_dh ! for CALL lim_thd_snwblow
44	#elif defined key_lim2
45	USE ice_2 , ONLY : u_ice, v_ice
46	USE par_ice_2 ! LIM-2 parameters
47	#endif
48	USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009)
49	USE sbcblk_algo_coare ! => turb_coare : COAREv3.0 (Fairall et al. 2003)
50	USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)
51	USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 31)
52	!
53	USE iom ! I/O manager library
54	USE in_out_manager ! I/O manager
55	USE lib_mpp ! distribued memory computing library
56	USE wrk_nemo ! work arrays
57	USE timing ! Timing
58	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
59	USE prtctl ! Print control
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_lim2 \|\| defined key_lim3
67	PUBLIC blk_ice_tau ! routine called in sbc_ice_lim module
68	PUBLIC blk_ice_flx ! routine called in sbc_ice_lim module
69	#endif
70
71	!!Lolo: should ultimately be moved in the module with all physical constants ?
72	!!gm : In principle, yes.
73	REAL(wp), PARAMETER :: Cp_dry = 1005.0 !: Specic heat of dry air, constant pressure [J/K/kg]
74	REAL(wp), PARAMETER :: Cp_vap = 1860.0 !: Specic heat of water vapor, constant pressure [J/K/kg]
75	REAL(wp), PARAMETER :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg]
76	REAL(wp), PARAMETER :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg]
77	REAL(wp), PARAMETER :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
78	REAL(wp), PARAMETER :: rctv0 = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
79
80	INTEGER , PARAMETER :: jpfld =10 ! maximum number of files to read
81	INTEGER , PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point
82	INTEGER , PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point
83	INTEGER , PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin)
84	INTEGER , PARAMETER :: jp_humi = 4 ! index of specific humidity ( % )
85	INTEGER , PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2)
86	INTEGER , PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2)
87	INTEGER , PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s)
88	INTEGER , PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s)
89	INTEGER , PARAMETER :: jp_slp = 9 ! index of snow (solid prcipitation) (kg/m2/s)
90	INTEGER , PARAMETER :: jp_tdif =10 ! index of tau diff associated to HF tau (N/m2) at T-point
91
92	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
93
94	! !!! Bulk parameters
95	REAL(wp), PARAMETER :: cpa = 1000.5 ! specific heat of air (only used for ice fluxes now...)
96	REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation
97	REAL(wp), PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant
98	REAL(wp), PARAMETER :: Cice = 1.4e-3 ! iovi 1.63e-3 ! transfer coefficient over ice
99	REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant
100	!
101	! !!* Namelist namsbc_blk : bulk parameters
102	LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008)
103	LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003)
104	LOGICAL :: ln_COARE_3p5 ! "COARE 3.5" algorithm (Edson et al. 2013)
105	LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 31)
106	!
107	LOGICAL :: ln_taudif ! logical flag to use the "mean of stress module - module of mean stress" data
108	REAL(wp) :: rn_pfac ! multiplication factor for precipitation
109	REAL(wp) :: rn_efac ! multiplication factor for evaporation (clem)
110	REAL(wp) :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
111	REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements
112	REAL(wp) :: rn_zu ! z(u) : height of wind measurements
113
114	INTEGER :: nblk ! choice of the bulk algorithm
115	! ! associated indices:
116	INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008)
117	INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003)
118	INTEGER, PARAMETER :: np_COARE_3p5 = 3 ! "COARE 3.5" algorithm (Edson et al. 2013)
119	INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 31)
120
121	!! * Substitutions
122	# include "vectopt_loop_substitute.h90"
123	!!----------------------------------------------------------------------
124	!! NEMO/OPA 3.7 , NEMO-consortium (2014)
125	!! $Id: sbcblk.F90 6416 2016-04-01 12:22:17Z clem $
126	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
127	!!----------------------------------------------------------------------
128	CONTAINS
129
130	SUBROUTINE sbc_blk_init
131	!!---------------------------------------------------------------------
132	!! * ROUTINE sbc_blk_init *
133	!!
134	!! ** Purpose : choose and initialize a bulk formulae formulation
135	!!
136	!! ** Method :
137	!!
138	!! C A U T I O N : never mask the surface stress fields
139	!! the stress is assumed to be in the (i,j) mesh referential
140	!!
141	!! ** Action :
142	!!
143	!!----------------------------------------------------------------------
144	INTEGER :: ifpr, jfld ! dummy loop indice and argument
145	INTEGER :: ios, ierror, ioptio ! Local integer
146	!!
147	CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files
148	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
149	TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read
150	TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " "
151	TYPE(FLD_N) :: sn_slp , sn_tdif ! " "
152	NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields
153	& sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif, &
154	& ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm
155	& cn_dir , ln_taudif, rn_zqt, rn_zu, rn_pfac, rn_efac, rn_vfac
156	!!---------------------------------------------------------------------
157	!
158	! !** read bulk namelist
159	REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters
160	READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
161	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist', lwp )
162	!
163	REWIND( numnam_cfg ) !* Namelist namsbc_blk in configuration namelist : bulk parameters
164	READ ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
165	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist', lwp )
166	!
167	IF(lwm) WRITE( numond, namsbc_blk )
168	!
169	! !** initialization of the chosen bulk formulae (+ check)
170	! !* select the bulk chosen in the namelist and check the choice
171	; ioptio = 0
172	IF( ln_NCAR ) THEN ; nblk = np_NCAR ; ioptio = ioptio + 1 ; ENDIF
173	IF( ln_COARE_3p0 ) THEN ; nblk = np_COARE_3p0 ; ioptio = ioptio + 1 ; ENDIF
174	IF( ln_COARE_3p5 ) THEN ; nblk = np_COARE_3p5 ; ioptio = ioptio + 1 ; ENDIF
175	IF( ln_ECMWF ) THEN ; nblk = np_ECMWF ; ioptio = ioptio + 1 ; ENDIF
176	!
177	IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
178	!
179	IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr
180	IF( sn_qsr%nfreqh /= 24 ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
181	IF( sn_qsr%ln_tint ) THEN
182	CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', &
183	& ' ==> We force time interpolation = .false. for qsr' )
184	sn_qsr%ln_tint = .false.
185	ENDIF
186	ENDIF
187	! !* set the bulk structure
188	! !- store namelist information in an array
189	slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj
190	slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw
191	slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi
192	slf_i(jp_prec) = sn_prec ; slf_i(jp_snow) = sn_snow
193	slf_i(jp_slp) = sn_slp ; slf_i(jp_tdif) = sn_tdif
194	!
195	lhftau = ln_taudif !- add an extra field if HF stress is used
196	jfld = jpfld - COUNT( (/.NOT.lhftau/) )
197	!
198	! !- allocate the bulk structure
199	ALLOCATE( sf(jfld), STAT=ierror )
200	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
201	DO ifpr= 1, jfld
202	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
203	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
204	END DO
205	! !- fill the bulk structure with namelist informations
206	CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
207	!
208	!
209	IF(lwp) THEN !** Control print
210	!
211	WRITE(numout,) ! namelist
212	WRITE(numout,*) ' Namelist namsbc_blk (other than data information):'
213	WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR
214	WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0
215	WRITE(numout,*) ' "COARE 3.5" algorithm (Edson et al. 2013) ln_COARE_3p5 = ', ln_COARE_3p0
216	WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 31) ln_ECMWF = ', ln_ECMWF
217	WRITE(numout,*) ' add High freq.contribution to the stress module ln_taudif = ', ln_taudif
218	WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt
219	WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu
220	WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac
221	WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac
222	WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac
223	WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))'
224	!
225	WRITE(numout,*)
226	SELECT CASE( nblk ) !* Print the choice of bulk algorithm
227	CASE( np_NCAR ) ; WRITE(numout,*) ' ===>> "NCAR" algorithm (Large and Yeager 2008)'
228	CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ===>> "COARE 3.0" algorithm (Fairall et al. 2003)'
229	CASE( np_COARE_3p5 ) ; WRITE(numout,*) ' ===>> "COARE 3.5" algorithm (Edson et al. 2013)'
230	CASE( np_ECMWF ) ; WRITE(numout,*) ' ===>> "ECMWF" algorithm (IFS cycle 31)'
231	END SELECT
232	!
233	ENDIF
234	!
235	END SUBROUTINE sbc_blk_init
236
237
238	SUBROUTINE sbc_blk( kt )
239	!!---------------------------------------------------------------------
240	!! * ROUTINE sbc_blk *
241	!!
242	!! ** Purpose : provide at each time step the surface ocean fluxes
243	!! (momentum, heat, freshwater and runoff)
244	!!
245	!! ** Method : (1) READ each fluxes in NetCDF files:
246	!! the 10m wind velocity (i-component) (m/s) at T-point
247	!! the 10m wind velocity (j-component) (m/s) at T-point
248	!! the 10m or 2m specific humidity ( % )
249	!! the solar heat (W/m2)
250	!! the Long wave (W/m2)
251	!! the 10m or 2m air temperature (Kelvin)
252	!! the total precipitation (rain+snow) (Kg/m2/s)
253	!! the snow (solid prcipitation) (kg/m2/s)
254	!! the tau diff associated to HF tau (N/m2) at T-point (ln_taudif=T)
255	!! (2) CALL blk_oce
256	!!
257	!! C A U T I O N : never mask the surface stress fields
258	!! the stress is assumed to be in the (i,j) mesh referential
259	!!
260	!! ** Action : defined at each time-step at the air-sea interface
261	!! - utau, vtau i- and j-component of the wind stress
262	!! - taum wind stress module at T-point
263	!! - wndm wind speed module at T-point over free ocean or leads in presence of sea-ice
264	!! - qns, qsr non-solar and solar heat fluxes
265	!! - emp upward mass flux (evapo. - precip.)
266	!! - sfx salt flux due to freezing/melting (non-zero only if ice is present)
267	!! (set in limsbc(_2).F90)
268	!!
269	!! ** References : Large & Yeager, 2004 / Large & Yeager, 2008
270	!! Brodeau et al. Ocean Modelling 2010
271	!!----------------------------------------------------------------------
272	INTEGER, INTENT(in) :: kt ! ocean time step
273	!!---------------------------------------------------------------------
274	!
275	CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step
276	!
277	! ! compute the surface ocean fluxes using bulk formulea
278	IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )
279
280	#if defined key_cice
281	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
282	qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1)
283	qsr_ice(:,:,1) = sf(jp_qsr )%fnow(:,:,1)
284	tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)
285	qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1)
286	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
287	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
288	wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
289	wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
290	ENDIF
291	#endif
292	!
293	END SUBROUTINE sbc_blk
294
295
296	SUBROUTINE blk_oce( kt, sf, pst, pu, pv )
297	!!---------------------------------------------------------------------
298	!! * ROUTINE blk_oce *
299	!!
300	!! ** Purpose : provide the momentum, heat and freshwater fluxes at
301	!! the ocean surface at each time step
302	!!
303	!! ** Method : bulk formulea for the ocean using atmospheric
304	!! fields read in sbc_read
305	!!
306	!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
307	!! - vtau : j-component of the stress at V-point (N/m2)
308	!! - taum : Wind stress module at T-point (N/m2)
309	!! - wndm : Wind speed module at T-point (m/s)
310	!! - qsr : Solar heat flux over the ocean (W/m2)
311	!! - qns : Non Solar heat flux over the ocean (W/m2)
312	!! - emp : evaporation minus precipitation (kg/m2/s)
313	!!
314	!! ** Nota : sf has to be a dummy argument for AGRIF on NEC
315	!!---------------------------------------------------------------------
316	INTEGER , INTENT(in ) :: kt ! time step index
317	TYPE(fld), INTENT(inout), DIMENSION(:) :: sf ! input data
318	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pst ! surface temperature [Celcius]
319	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s]
320	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s]
321	!
322	INTEGER :: ji, jj ! dummy loop indices
323	REAL(wp) :: zztmp ! local variable
324	REAL(wp), DIMENSION(:,:), POINTER :: zwnd_i, zwnd_j ! wind speed components at T-point
325	REAL(wp), DIMENSION(:,:), POINTER :: zsq ! specific humidity at pst
326	REAL(wp), DIMENSION(:,:), POINTER :: zqlw, zqsb ! long wave and sensible heat fluxes
327	REAL(wp), DIMENSION(:,:), POINTER :: zqla, zevap ! latent heat fluxes and evaporation
328	REAL(wp), DIMENSION(:,:), POINTER :: Cd ! transfer coefficient for momentum (tau)
329	REAL(wp), DIMENSION(:,:), POINTER :: Ch ! transfer coefficient for sensible heat (Q_sens)
330	REAL(wp), DIMENSION(:,:), POINTER :: Ce ! tansfert coefficient for evaporation (Q_lat)
331	REAL(wp), DIMENSION(:,:), POINTER :: zst ! surface temperature in Kelvin
332	REAL(wp), DIMENSION(:,:), POINTER :: zt_zu ! air temperature at wind speed height
333	REAL(wp), DIMENSION(:,:), POINTER :: zq_zu ! air spec. hum. at wind speed height
334	REAL(wp), DIMENSION(:,:), POINTER :: zU_zu ! bulk wind speed at height zu [m/s]
335	REAL(wp), DIMENSION(:,:), POINTER :: ztpot ! potential temperature of air at z=rn_zqt [K]
336	REAL(wp), DIMENSION(:,:), POINTER :: zrhoa ! density of air [kg/m^3]
337	!!---------------------------------------------------------------------
338	!
339	IF( nn_timing == 1 ) CALL timing_start('blk_oce')
340	!
341	CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap )
342	CALL wrk_alloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )
343	CALL wrk_alloc( jpi,jpj, zU_zu, ztpot, zrhoa )
344	!
345
346	! local scalars ( place there for vector optimisation purposes)
347	zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
348
349	! ----------------------------------------------------------------------------- !
350	! 0 Wind components and module at T-point relative to the moving ocean !
351	! ----------------------------------------------------------------------------- !
352
353	! ... components ( U10m - U_oce ) at T-point (unmasked)
354	!!gm move zwnd_i (_j) set to zero inside the key_cyclone ???
355	zwnd_i(:,:) = 0._wp
356	zwnd_j(:,:) = 0._wp
357	#if defined key_cyclone
358	CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012)
359	DO jj = 2, jpjm1
360	DO ji = fs_2, fs_jpim1 ! vect. opt.
361	sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
362	sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
363	END DO
364	END DO
365	#endif
366	DO jj = 2, jpjm1
367	DO ji = fs_2, fs_jpim1 ! vect. opt.
368	zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) )
369	zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) )
370	END DO
371	END DO
372	CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. )
373	CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. )
374	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
375	wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) &
376	& + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1)
377
378	! ----------------------------------------------------------------------------- !
379	! I Radiative FLUXES !
380	! ----------------------------------------------------------------------------- !
381
382	! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
383	zztmp = 1. - albo
384	IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
385	ELSE ; qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1)
386	ENDIF
387
388	zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)zst(:,:)zst(:,:)zst(:,:) ) tmask(:,:,1) ! Long Wave
389
390
391
392	! ----------------------------------------------------------------------------- !
393	! II Turbulent FLUXES !
394	! ----------------------------------------------------------------------------- !
395
396	! ... specific humidity at SST and IST tmask(
397	zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
398	!!
399	!! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
400	!! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
401	!! (since reanalysis products provide T at z, not theta !)
402	ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt
403
404	SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point
405	!
406	CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2
407	& Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
408	CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0
409	& Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
410	CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5
411	& Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
412	CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF
413	& Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
414	CASE DEFAULT
415	CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
416	END SELECT
417
418	! ! Compute true air density :
419	IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrhogravrn_zu from SLP...)
420	zrhoa(:,:) = rho_air( zt_zu(:,:) , zq_zu(:,:) , sf(jp_slp)%fnow(:,:,1) )
421	ELSE ! At zt:
422	zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
423	END IF
424
425	DO jj = 1, jpj ! tau module, i and j component
426	DO ji = 1, jpi
427	zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd(ji,jj) ! using bulk wind speed
428	taum (ji,jj) = zztmp * wndm (ji,jj)
429	zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
430	zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
431	END DO
432	END DO
433
434	! ! add the HF tau contribution to the wind stress module
435	IF( lhftau ) taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
436
437	CALL iom_put( "taum_oce", taum ) ! output wind stress module
438
439	! ... utau, vtau at U- and V_points, resp.
440	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
441	! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
442	DO jj = 1, jpjm1
443	DO ji = 1, fs_jpim1
444	utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) &
445	& * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
446	vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) &
447	& * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
448	END DO
449	END DO
450	CALL lbc_lnk( utau(:,:), 'U', -1. )
451	CALL lbc_lnk( vtau(:,:), 'V', -1. )
452
453
454	! Turbulent fluxes over ocean
455	! -----------------------------
456
457	! zqla used as temporary array, for rho*U (common term of bulk formulae):
458	zqla(:,:) = zrhoa(:,:) * zU_zu(:,:)
459
460	IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
461	!! q_air and t_air are given at 10m (wind reference height)
462	zevap(:,:) = rn_efacMAX( 0._wp, zqla(:,:)Ce(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed
463	zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))zqla(:,:)Ch(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed
464	ELSE
465	!! q_air and t_air are not given at 10m (wind reference height)
466	! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
467	zevap(:,:) = rn_efacMAX( 0._wp, zqla(:,:)Ce(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed
468	zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))zqla(:,:)Ch(:,:)*(zst(:,:) - zt_zu(:,:) ) ! Sensible Heat ! using bulk wind speed
469	ENDIF
470
471	zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:) ! Latent Heat flux
472
473
474	IF(ln_ctl) THEN
475	CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce , clinfo2=' Ce : ' )
476	CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=Ch , clinfo2=' Ch : ' )
477	CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' )
478	CALL prt_ctl( tab2d_1=zsq , clinfo1=' blk_oce: zsq : ', tab2d_2=zst, clinfo2=' zst : ' )
479	CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, &
480	& tab2d_2=vtau , clinfo2= ' vtau : ', mask2=vmask )
481	CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce: wndm : ')
482	CALL prt_ctl( tab2d_1=zst , clinfo1=' blk_oce: zst : ')
483	ENDIF
484
485	! ----------------------------------------------------------------------------- !
486	! III Total FLUXES !
487	! ----------------------------------------------------------------------------- !
488	!
489	emp (:,:) = ( zevap(:,:) & ! mass flux (evap. - precip.)
490	& - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1)
491	!
492	qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar
493	& - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus & ! remove latent melting heat for solid precip
494	& - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST
495	& + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac & ! add liquid precip heat content at Tair
496	& * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp &
497	& + sf(jp_snow)%fnow(:,:,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)
498	& * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1)
499	!
500	#if defined key_lim3
501	qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) ! non solar without emp (only needed by LIM3)
502	qsr_oce(:,:) = qsr(:,:)
503	#endif
504	!
505	IF ( nn_ice == 0 ) THEN
506	CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean
507	CALL iom_put( "qsb_oce" , - zqsb ) ! output downward sensible heat over the ocean
508	CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean
509	CALL iom_put( "qemp_oce", qns-zqlw+zqsb+zqla ) ! output downward heat content of E-P over the ocean
510	CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean
511	CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean
512	CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean
513	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac ! output total precipitation [kg/m2/s]
514	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac ! output solid precipitation [kg/m2/s]
515	CALL iom_put( 'snowpre', sprecip * 86400. ) ! Snow
516	CALL iom_put( 'precip' , tprecip * 86400. ) ! Total precipitation
517	ENDIF
518	!
519	IF(ln_ctl) THEN
520	CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=zqlw , clinfo2=' zqlw : ')
521	CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ')
522	CALL prt_ctl(tab2d_1=pst , clinfo1=' blk_oce: pst : ', tab2d_2=emp , clinfo2=' emp : ')
523	CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, &
524	& tab2d_2=vtau , clinfo2= ' vtau : ' , mask2=vmask )
525	ENDIF
526	!
527	CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap )
528	CALL wrk_dealloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )
529	CALL wrk_dealloc( jpi,jpj, zU_zu, ztpot, zrhoa )
530	!
531	IF( nn_timing == 1 ) CALL timing_stop('blk_oce')
532	!
533	END SUBROUTINE blk_oce
534
535	#if defined key_lim2 \|\| defined key_lim3
536
537	SUBROUTINE blk_ice_tau
538	!!---------------------------------------------------------------------
539	!! * ROUTINE blk_ice_tau *
540	!!
541	!! ** Purpose : provide the surface boundary condition over sea-ice
542	!!
543	!! ** Method : compute momentum using bulk formulation
544	!! formulea, ice variables and read atmospheric fields.
545	!! NB: ice drag coefficient is assumed to be a constant
546	!!---------------------------------------------------------------------
547	INTEGER :: ji, jj ! dummy loop indices
548	!
549	REAL(wp), DIMENSION(:,:) , POINTER :: zcoef_wnorm
550	!
551	REAL(wp) :: zwnorm_f, zwndi_f , zwndj_f ! relative wind module and components at F-point
552	REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point
553	!!---------------------------------------------------------------------
554	!
555	IF( nn_timing == 1 ) CALL timing_start('blk_ice_tau')
556	!
557	CALL wrk_alloc( jpi,jpj, zcoef_wnorm )
558
559	! local scalars ( place there for vector optimisation purposes)
560	! Computing density of air! Way denser that 1.2 over sea-ice !!!
561	!!
562	zcoef_wnorm (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
563	zcoef_wnorm (:,:) = Cice * zcoef_wnorm (:,:)
564
565	!!gm brutal....
566	utau_ice (:,:) = 0._wp
567	vtau_ice (:,:) = 0._wp
568	wndm_ice (:,:) = 0._wp
569	!!gm end
570
571	! ----------------------------------------------------------------------------- !
572	! Wind components and module relative to the moving ocean ( U10m - U_ice ) !
573	! ----------------------------------------------------------------------------- !
574	SELECT CASE( cp_ice_msh )
575	CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation)
576	! and scalar wind at T-point ( = \| U10m - U_ice \| ) (masked)
577	DO jj = 2, jpjm1
578	DO ji = 2, jpim1 ! B grid : NO vector opt
579	! ... scalar wind at I-point (fld being at T-point)
580	zwndi_f = 0.25 * ( sf(jp_wndi)%fnow(ji-1,jj ,1) + sf(jp_wndi)%fnow(ji ,jj ,1) &
581	& + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji ,jj-1,1) ) - rn_vfac * u_ice(ji,jj)
582	zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) &
583	& + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * v_ice(ji,jj)
584	zwnorm_f = zcoef_wnorm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
585	! ... ice stress at I-point
586	utau_ice(ji,jj) = zwnorm_f * zwndi_f
587	vtau_ice(ji,jj) = zwnorm_f * zwndj_f
588	! ... scalar wind at T-point (fld being at T-point)
589	zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( u_ice(ji,jj+1) + u_ice(ji+1,jj+1) &
590	& + u_ice(ji,jj ) + u_ice(ji+1,jj ) )
591	zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) &
592	& + v_ice(ji,jj ) + v_ice(ji+1,jj ) )
593	wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
594	END DO
595	END DO
596	CALL lbc_lnk( utau_ice, 'I', -1. )
597	CALL lbc_lnk( vtau_ice, 'I', -1. )
598	CALL lbc_lnk( wndm_ice, 'T', 1. )
599	!
600	CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean)
601	DO jj = 2, jpj
602	DO ji = fs_2, jpi ! vect. opt.
603	zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) )
604	zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )
605	wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
606	END DO
607	END DO
608	DO jj = 2, jpjm1
609	DO ji = fs_2, fs_jpim1 ! vect. opt.
610	utau_ice(ji,jj) = 0.5zcoef_wnorm(ji,jj) ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) &
611	& * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
612	vtau_ice(ji,jj) = 0.5zcoef_wnorm(ji,jj) ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) &
613	& * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
614	END DO
615	END DO
616	CALL lbc_lnk( utau_ice, 'U', -1. )
617	CALL lbc_lnk( vtau_ice, 'V', -1. )
618	CALL lbc_lnk( wndm_ice, 'T', 1. )
619	!
620	END SELECT
621
622	IF(ln_ctl) THEN
623	CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ')
624	CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
625	ENDIF
626
627	CALL wrk_dealloc( jpi,jpj, zcoef_wnorm )
628
629	IF( nn_timing == 1 ) CALL timing_stop('blk_ice_tau')
630
631	END SUBROUTINE blk_ice_tau
632
633
634	SUBROUTINE blk_ice_flx( ptsu, palb )
635	!!---------------------------------------------------------------------
636	!! * ROUTINE blk_ice_flx *
637	!!
638	!! ** Purpose : provide the surface boundary condition over sea-ice
639	!!
640	!! ** Method : compute heat and freshwater exchanged
641	!! between atmosphere and sea-ice using bulk formulation
642	!! formulea, ice variables and read atmmospheric fields.
643	!!
644	!! caution : the net upward water flux has with mm/day unit
645	!!---------------------------------------------------------------------
646	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature
647	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies)
648	!!
649	INTEGER :: ji, jj, jl ! dummy loop indices
650	REAL(wp) :: zst2, zst3 ! local variable
651	REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
652	REAL(wp) :: zztmp, z1_lsub ! - -
653	REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice
654	REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice
655	REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqlw ! long wave heat sensitivity over ice
656	REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqsb ! sensible heat sensitivity over ice
657	REAL(wp), DIMENSION(:,:) , POINTER :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3)
658	REAL(wp), DIMENSION(:,:) , POINTER :: zrhoa
659	!!---------------------------------------------------------------------
660	!
661	IF( nn_timing == 1 ) CALL timing_start('blk_ice_flx')
662	!
663	CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )
664	CALL wrk_alloc( jpi,jpj, zrhoa)
665	!
666	! local scalars ( place there for vector optimisation purposes)
667	zcoef_dqlw = 4.0 * 0.95 * Stef
668	zcoef_dqla = -Ls * Cice * 11637800. * (-5897.8)
669	!
670	zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
671	!
672	zztmp = 1. / ( 1. - albo )
673	! ! ========================== !
674	DO jl = 1, jpl ! Loop over ice categories !
675	! ! ========================== !
676	DO jj = 1 , jpj
677	DO ji = 1, jpi
678	! ----------------------------!
679	! I Radiative FLUXES !
680	! ----------------------------!
681	zst2 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
682	zst3 = ptsu(ji,jj,jl) * zst2
683	! Short Wave (sw)
684	qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
685	! Long Wave (lw)
686	z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
687	! lw sensitivity
688	z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
689
690	! ----------------------------!
691	! II Turbulent FLUXES !
692	! ----------------------------!
693
694	! ... turbulent heat fluxes
695	! Sensible Heat
696	z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
697	! Latent Heat
698	qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Cice * wndm_ice(ji,jj) &
699	& * ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
700	! Latent heat sensitivity for ice (Dqla/Dt)
701	IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
702	dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) )
703	ELSE
704	dqla_ice(ji,jj,jl) = 0._wp
705	ENDIF
706
707	! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
708	z_dqsb(ji,jj,jl) = zrhoa(ji,jj)cpaCice * wndm_ice(ji,jj)
709
710	! ----------------------------!
711	! III Total FLUXES !
712	! ----------------------------!
713	! Downward Non Solar flux
714	qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
715	! Total non solar heat flux sensitivity for ice
716	dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
717	END DO
718	!
719	END DO
720	!
721	END DO
722	!
723	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac ! total precipitation [kg/m2/s]
724	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac ! solid precipitation [kg/m2/s]
725	CALL iom_put( 'snowpre', sprecip * 86400. ) ! Snow precipitation
726	CALL iom_put( 'precip' , tprecip * 86400. ) ! Total precipitation
727
728	#if defined key_lim3
729	CALL wrk_alloc( jpi,jpj, zevap, zsnw )
730
731	! --- evaporation --- !
732	z1_lsub = 1._wp / Lsub
733	evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub ! sublimation
734	devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub ! d(sublimation)/dT
735	zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean
736
737	! --- evaporation minus precipitation --- !
738	zsnw(:,:) = 0._wp
739	CALL lim_thd_snwblow( pfrld, zsnw ) ! snow distribution over ice after wind blowing
740	emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
741	emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
742	emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
743
744	! --- heat flux associated with emp --- !
745	qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst
746	& + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & ! liquid precip at Tair
747	& + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow)
748	& ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
749	qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only)
750	& ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
751
752	! --- total solar and non solar fluxes --- !
753	qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:)
754	qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
755
756	! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
757	qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
758
759	CALL wrk_dealloc( jpi,jpj, zevap, zsnw )
760	#endif
761
762	!--------------------------------------------------------------------
763	! FRACTIONs of net shortwave radiation which is not absorbed in the
764	! thin surface layer and penetrates inside the ice cover
765	! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
766	!
767	fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
768	fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
769	!
770	!
771	IF(ln_ctl) THEN
772	CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl)
773	CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
774	CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl)
775	CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl)
776	CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl)
777	CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ')
778	ENDIF
779
780	CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )
781	CALL wrk_dealloc( jpi,jpj, zrhoa )
782	!
783	IF( nn_timing == 1 ) CALL timing_stop('blk_ice_flx')
784
785	END SUBROUTINE blk_ice_flx
786
787	#endif
788
789	FUNCTION rho_air( ptak, pqa, pslp )
790	!!-------------------------------------------------------------------------------
791	!! * FUNCTION rho_air *
792	!!
793	!! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
794	!!
795	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
796	!!-------------------------------------------------------------------------------
797	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
798	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
799	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa]
800	REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3]
801	!!-------------------------------------------------------------------------------
802	!
803	rho_air = pslp / ( R_dryptak ( 1._wp + rctv0*pqa ) )
804	!
805	END FUNCTION rho_air
806
807
808	FUNCTION cp_air( pqa )
809	!!-------------------------------------------------------------------------------
810	!! * FUNCTION cp_air *
811	!!
812	!! ** Purpose : Compute specific heat (Cp) of moist air
813	!!
814	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
815	!!-------------------------------------------------------------------------------
816	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
817	REAL(wp), DIMENSION(jpi,jpj) :: cp_air ! specific heat of moist air [J/K/kg]
818	!!-------------------------------------------------------------------------------
819	!
820	Cp_air = Cp_dry + Cp_vap * pqa
821	!
822	END FUNCTION cp_air
823
824
825	FUNCTION q_sat( ptak, pslp )
826	!!----------------------------------------------------------------------------------
827	!! * FUNCTION q_sat *
828	!!
829	!! ** Purpose : Specific humidity at saturation in [kg/kg]
830	!! Based on accurate estimate of "e_sat"
831	!! aka saturation water vapor (Goff, 1957)
832	!!
833	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
834	!!----------------------------------------------------------------------------------
835	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
836	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa]
837	REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg]
838	!
839	INTEGER :: ji, jj ! dummy loop indices
840	REAL(wp) :: ze_sat, ztmp ! local scalar
841	!!----------------------------------------------------------------------------------
842	!
843	DO jj = 1, jpj
844	DO ji = 1, jpi
845	!
846	ztmp = rt0 / ptak(ji,jj)
847	!
848	! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
849	ze_sat = 10.*( 10.79574(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) &
850	& + 1.5047510.(-4)(1. - 10.*(-8.2969(ptak(ji,jj)/rt0 - 1.)) ) &
851	& + 0.4287310.(-3)(10.*(4.76955(1. - ztmp)) - 1.) + 0.78614 )
852	!
853	q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wppslp(ji,jj) - (1._wp - reps0)ze_sat ) ! 0.01 because SLP is in [Pa]
854	!
855	END DO
856	END DO
857	!
858	END FUNCTION q_sat
859
860
861	FUNCTION gamma_moist( ptak, pqa )
862	!!----------------------------------------------------------------------------------
863	!! * FUNCTION gamma_moist *
864	!!
865	!! ** Purpose : Compute the moist adiabatic lapse-rate.
866	!! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
867	!! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
868	!!
869	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
870	!!----------------------------------------------------------------------------------
871	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
872	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg]
873	REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate
874	!
875	INTEGER :: ji, jj ! dummy loop indices
876	REAL(wp) :: zrv, ziRT ! local scalar
877	!!----------------------------------------------------------------------------------
878	!
879	DO jj = 1, jpj
880	DO ji = 1, jpi
881	zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
882	ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT
883	gamma_moist(ji,jj) = grav * ( 1. + cevapzrvziRT ) / ( Cp_dry + cevapcevapzrvreps0ziRT/ptak(ji,jj) )
884	END DO
885	END DO
886	!
887	END FUNCTION gamma_moist
888
889
890	FUNCTION L_vap( psst )
891	!!---------------------------------------------------------------------------------
892	!! * FUNCTION L_vap *
893	!!
894	!! ** Purpose : Compute the latent heat of vaporization of water from temperature
895	!!
896	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
897	!!----------------------------------------------------------------------------------
898	REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg]
899	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]
900	!!----------------------------------------------------------------------------------
901	!
902	L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6
903	!
904	END FUNCTION L_vap
905
906	!!======================================================================
907	END MODULE sbcblk

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

Download in other formats: