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_core.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 10 years ago
remove remaining bugs in LIM3, so that it can run in both regional and global config
Property svn:keywords set to `Id`
File size: 58.4 KB

Line
1	MODULE sbcblk_core
2	!!======================================================================
3	!! * MODULE sbcblk_core *
4	!! Ocean forcing: momentum, heat and freshwater flux formulation
5	!!=====================================================================
6	!! History : 1.0 ! 2004-08 (U. Schweckendiek) Original code
7	!! 2.0 ! 2005-04 (L. Brodeau, A.M. Treguier) additions:
8	!! - new bulk routine for efficiency
9	!! - WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files !!!!
10	!! - file names and file characteristics in namelist
11	!! - Implement reading of 6-hourly fields
12	!! 3.0 ! 2006-06 (G. Madec) sbc rewritting
13	!! - ! 2006-12 (L. Brodeau) Original code for TURB_CORE_2Z
14	!! 3.2 ! 2009-04 (B. Lemaire) Introduce iom_put
15	!! 3.3 ! 2010-10 (S. Masson) add diurnal cycle
16	!! 3.4 ! 2011-11 (C. Harris) Fill arrays required by CICE
17	!!----------------------------------------------------------------------
18
19	!!----------------------------------------------------------------------
20	!! sbc_blk_core : bulk formulation as ocean surface boundary condition
21	!! (forced mode, CORE bulk formulea)
22	!! blk_oce_core : ocean: computes momentum, heat and freshwater fluxes
23	!! blk_ice_core : ice : computes momentum, heat and freshwater fluxes
24	!! turb_core : computes the CORE turbulent transfer coefficients
25	!!----------------------------------------------------------------------
26	USE oce ! ocean dynamics and tracers
27	USE dom_oce ! ocean space and time domain
28	USE phycst ! physical constants
29	USE fldread ! read input fields
30	USE sbc_oce ! Surface boundary condition: ocean fields
31	USE cyclone ! Cyclone 10m wind form trac of cyclone centres
32	USE sbcdcy ! surface boundary condition: diurnal cycle
33	USE iom ! I/O manager library
34	USE in_out_manager ! I/O manager
35	USE lib_mpp ! distribued memory computing library
36	USE wrk_nemo ! work arrays
37	USE timing ! Timing
38	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
39	USE prtctl ! Print control
40	USE sbcwave,ONLY : cdn_wave !wave module
41	#if defined key_lim3 \|\| defined key_cice
42	USE sbc_ice ! Surface boundary condition: ice fields
43	#endif
44	USE lib_fortran ! to use key_nosignedzero
45
46	IMPLICIT NONE
47	PRIVATE
48
49	PUBLIC sbc_blk_core ! routine called in sbcmod module
50	PUBLIC blk_ice_core ! routine called in sbc_ice_lim module
51	PUBLIC blk_ice_meanqsr ! routine called in sbc_ice_lim module
52	PUBLIC turb_core_2z ! routine calles in sbcblk_mfs module
53
54	INTEGER , PARAMETER :: jpfld = 9 ! maximum number of files to read
55	INTEGER , PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point
56	INTEGER , PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point
57	INTEGER , PARAMETER :: jp_humi = 3 ! index of specific humidity ( % )
58	INTEGER , PARAMETER :: jp_qsr = 4 ! index of solar heat (W/m2)
59	INTEGER , PARAMETER :: jp_qlw = 5 ! index of Long wave (W/m2)
60	INTEGER , PARAMETER :: jp_tair = 6 ! index of 10m air temperature (Kelvin)
61	INTEGER , PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s)
62	INTEGER , PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s)
63	INTEGER , PARAMETER :: jp_tdif = 9 ! index of tau diff associated to HF tau (N/m2) at T-point
64
65	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
66
67	! !!! CORE bulk parameters
68	REAL(wp), PARAMETER :: rhoa = 1.22 ! air density
69	REAL(wp), PARAMETER :: cpa = 1000.5 ! specific heat of air
70	REAL(wp), PARAMETER :: Lv = 2.5e6 ! latent heat of vaporization
71	REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation
72	REAL(wp), PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant
73	REAL(wp), PARAMETER :: Cice = 1.4e-3 ! iovi 1.63e-3 ! transfer coefficient over ice
74	REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant
75
76	! !!* Namelist namsbc_core : CORE bulk parameters
77	LOGICAL :: ln_2m ! logical flag for height of air temp. and hum
78	LOGICAL :: ln_taudif ! logical flag to use the "mean of stress module - module of mean stress" data
79	REAL(wp) :: rn_pfac ! multiplication factor for precipitation
80	REAL(wp) :: rn_efac ! multiplication factor for evaporation (clem)
81	REAL(wp) :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
82	LOGICAL :: ln_bulk2z ! logical flag for case where z(q,t) and z(u) are specified in the namelist
83	REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements
84	REAL(wp) :: rn_zu ! z(u) : height of wind measurements
85
86	!! * Substitutions
87	# include "domzgr_substitute.h90"
88	# include "vectopt_loop_substitute.h90"
89	!!----------------------------------------------------------------------
90	!! NEMO/OPA 3.3 , NEMO-consortium (2010)
91	!! $Id$
92	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
93	!!----------------------------------------------------------------------
94	CONTAINS
95
96	SUBROUTINE sbc_blk_core( kt )
97	!!---------------------------------------------------------------------
98	!! * ROUTINE sbc_blk_core *
99	!!
100	!! ** Purpose : provide at each time step the surface ocean fluxes
101	!! (momentum, heat, freshwater and runoff)
102	!!
103	!! ** Method : (1) READ each fluxes in NetCDF files:
104	!! the 10m wind velocity (i-component) (m/s) at T-point
105	!! the 10m wind velocity (j-component) (m/s) at T-point
106	!! the 10m or 2m specific humidity ( % )
107	!! the solar heat (W/m2)
108	!! the Long wave (W/m2)
109	!! the 10m or 2m air temperature (Kelvin)
110	!! the total precipitation (rain+snow) (Kg/m2/s)
111	!! the snow (solid prcipitation) (kg/m2/s)
112	!! the tau diff associated to HF tau (N/m2) at T-point (ln_taudif=T)
113	!! (2) CALL blk_oce_core
114	!!
115	!! C A U T I O N : never mask the surface stress fields
116	!! the stress is assumed to be in the (i,j) mesh referential
117	!!
118	!! ** Action : defined at each time-step at the air-sea interface
119	!! - utau, vtau i- and j-component of the wind stress
120	!! - taum, wndm wind stress and 10m wind modules at T-point
121	!! - qns, qsr non-solar and solar heat fluxes
122	!! - emp upward mass flux (evapo. - precip.)
123	!! - sfx salt flux due to freezing/melting (non-zero only if ice is present)
124	!! (set in limsbc(_2).F90)
125	!!----------------------------------------------------------------------
126	INTEGER, INTENT(in) :: kt ! ocean time step
127	!!
128	INTEGER :: ierror ! return error code
129	INTEGER :: ifpr ! dummy loop indice
130	INTEGER :: jfld ! dummy loop arguments
131	INTEGER :: ios ! Local integer output status for namelist read
132	!!
133	CHARACTER(len=100) :: cn_dir ! Root directory for location of core files
134	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
135	TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read
136	TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " "
137	TYPE(FLD_N) :: sn_tdif ! " "
138	NAMELIST/namsbc_core/ cn_dir , ln_2m , ln_taudif, rn_pfac, rn_efac, rn_vfac, &
139	& sn_wndi, sn_wndj, sn_humi , sn_qsr , &
140	& sn_qlw , sn_tair, sn_prec , sn_snow, &
141	& sn_tdif, rn_zqt , ln_bulk2z, rn_zu
142	!!---------------------------------------------------------------------
143
144	! ! ====================== !
145	IF( kt == nit000 ) THEN ! First call kt=nit000 !
146	! ! ====================== !
147	!
148	REWIND( numnam_ref ) ! Namelist namsbc_core in reference namelist : CORE bulk parameters
149	READ ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901)
150	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_core in reference namelist', lwp )
151
152	REWIND( numnam_cfg ) ! Namelist namsbc_core in configuration namelist : CORE bulk parameters
153	READ ( numnam_cfg, namsbc_core, IOSTAT = ios, ERR = 902 )
154	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_core in configuration namelist', lwp )
155
156	WRITE ( numond, namsbc_core )
157	! ! check: do we plan to use ln_dm2dc with non-daily forcing?
158	IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 ) &
159	& CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' )
160	IF( ln_dm2dc .AND. sn_qsr%ln_tint ) THEN
161	CALL ctl_warn( 'sbc_blk_core: ln_dm2dc is taking care of the temporal interpolation of daily qsr', &
162	& ' ==> We force time interpolation = .false. for qsr' )
163	sn_qsr%ln_tint = .false.
164	ENDIF
165	! ! store namelist information in an array
166	slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj
167	slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw
168	slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi
169	slf_i(jp_prec) = sn_prec ; slf_i(jp_snow) = sn_snow
170	slf_i(jp_tdif) = sn_tdif
171	!
172	lhftau = ln_taudif ! do we use HF tau information?
173	jfld = jpfld - COUNT( (/.NOT. lhftau/) )
174	!
175	ALLOCATE( sf(jfld), STAT=ierror ) ! set sf structure
176	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_core: unable to allocate sf structure' )
177	DO ifpr= 1, jfld
178	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
179	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
180	END DO
181	! ! fill sf with slf_i and control print
182	CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' )
183	!
184	sfx(:,:) = 0._wp ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
185	!
186	ENDIF
187
188	CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step
189
190	! ! compute the surface ocean fluxes using CORE bulk formulea
191	IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce_core( kt, sf, sst_m, ssu_m, ssv_m )
192
193	! If diurnal cycle is activated, compute a daily mean short waves flux for biogeochemistery
194	IF( ltrcdm2dc ) CALL blk_bio_meanqsr
195
196	#if defined key_cice
197	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
198	qlw_ice(:,:,1) = sf(jp_qlw)%fnow(:,:,1)
199	qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1)
200	tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)
201	qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1)
202	tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
203	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
204	wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
205	wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
206	ENDIF
207	#endif
208	!
209	END SUBROUTINE sbc_blk_core
210
211
212	SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv )
213	!!---------------------------------------------------------------------
214	!! * ROUTINE blk_core *
215	!!
216	!! ** Purpose : provide the momentum, heat and freshwater fluxes at
217	!! the ocean surface at each time step
218	!!
219	!! ** Method : CORE bulk formulea for the ocean using atmospheric
220	!! fields read in sbc_read
221	!!
222	!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
223	!! - vtau : j-component of the stress at V-point (N/m2)
224	!! - taum : Wind stress module at T-point (N/m2)
225	!! - wndm : Wind speed module at T-point (m/s)
226	!! - qsr : Solar heat flux over the ocean (W/m2)
227	!! - qns : Non Solar heat flux over the ocean (W/m2)
228	!! - evap : Evaporation over the ocean (kg/m2/s)
229	!! - emp : evaporation minus precipitation (kg/m2/s)
230	!!
231	!! ** Nota : sf has to be a dummy argument for AGRIF on NEC
232	!!---------------------------------------------------------------------
233	INTEGER , INTENT(in ) :: kt ! time step index
234	TYPE(fld), INTENT(inout), DIMENSION(:) :: sf ! input data
235	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pst ! surface temperature [Celcius]
236	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s]
237	REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s]
238	!
239	INTEGER :: ji, jj ! dummy loop indices
240	REAL(wp) :: zcoef_qsatw, zztmp ! local variable
241	REAL(wp), DIMENSION(:,:), POINTER :: zwnd_i, zwnd_j ! wind speed components at T-point
242	REAL(wp), DIMENSION(:,:), POINTER :: zqsatw ! specific humidity at pst
243	REAL(wp), DIMENSION(:,:), POINTER :: zqlw, zqsb ! long wave and sensible heat fluxes
244	REAL(wp), DIMENSION(:,:), POINTER :: zqla, zevap ! latent heat fluxes and evaporation
245	REAL(wp), DIMENSION(:,:), POINTER :: Cd ! transfer coefficient for momentum (tau)
246	REAL(wp), DIMENSION(:,:), POINTER :: Ch ! transfer coefficient for sensible heat (Q_sens)
247	REAL(wp), DIMENSION(:,:), POINTER :: Ce ! tansfert coefficient for evaporation (Q_lat)
248	REAL(wp), DIMENSION(:,:), POINTER :: zst ! surface temperature in Kelvin
249	REAL(wp), DIMENSION(:,:), POINTER :: zt_zu ! air temperature at wind speed height
250	REAL(wp), DIMENSION(:,:), POINTER :: zq_zu ! air spec. hum. at wind speed height
251	!!---------------------------------------------------------------------
252	!
253	IF( nn_timing == 1 ) CALL timing_start('blk_oce_core')
254	!
255	CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap )
256	CALL wrk_alloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )
257	!
258	! local scalars ( place there for vector optimisation purposes)
259	zcoef_qsatw = 0.98 * 640380. / rhoa
260
261	zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
262
263	! ----------------------------------------------------------------------------- !
264	! 0 Wind components and module at T-point relative to the moving ocean !
265	! ----------------------------------------------------------------------------- !
266
267	! ... components ( U10m - U_oce ) at T-point (unmasked)
268	zwnd_i(:,:) = 0.e0
269	zwnd_j(:,:) = 0.e0
270	#if defined key_cyclone
271	# if defined key_vectopt_loop
272	!CDIR COLLAPSE
273	# endif
274	CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add Manu !
275	DO jj = 2, jpjm1
276	DO ji = fs_2, fs_jpim1 ! vect. opt.
277	sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
278	sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
279	END DO
280	END DO
281	#endif
282	#if defined key_vectopt_loop
283	!CDIR COLLAPSE
284	#endif
285	DO jj = 2, jpjm1
286	DO ji = fs_2, fs_jpim1 ! vect. opt.
287	zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) )
288	zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) )
289	END DO
290	END DO
291	CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. )
292	CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. )
293	! ... scalar wind ( = \| U10m - U_oce \| ) at T-point (masked)
294	!CDIR NOVERRCHK
295	!CDIR COLLAPSE
296	wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) &
297	& + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1)
298
299	! ----------------------------------------------------------------------------- !
300	! I Radiative FLUXES !
301	! ----------------------------------------------------------------------------- !
302
303	! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
304	zztmp = 1. - albo
305	IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
306	ELSE ; qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1)
307	ENDIF
308	!CDIR COLLAPSE
309	zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)zst(:,:)zst(:,:)zst(:,:) ) tmask(:,:,1) ! Long Wave
310	! ----------------------------------------------------------------------------- !
311	! II Turbulent FLUXES !
312	! ----------------------------------------------------------------------------- !
313
314	! ... specific humidity at SST and IST
315	!CDIR NOVERRCHK
316	!CDIR COLLAPSE
317	zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) )
318
319	! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point :
320	IF( ln_2m ) THEN
321	!! If air temp. and spec. hum. are given at different height (2m) than wind (10m) :
322	CALL TURB_CORE_2Z(2.,10., zst , sf(jp_tair)%fnow, &
323	& zqsatw, sf(jp_humi)%fnow, wndm, &
324	& Cd , Ch , Ce , &
325	& zt_zu , zq_zu )
326	ELSE IF( ln_bulk2z ) THEN
327	!! If the height of the air temp./spec. hum. and wind are to be specified by hand :
328	IF( rn_zqt == rn_zu ) THEN
329	!! If air temp. and spec. hum. are at the same height as wind :
330	CALL TURB_CORE_1Z( rn_zu, zst , sf(jp_tair)%fnow(:,:,1), &
331	& zqsatw, sf(jp_humi)%fnow(:,:,1), wndm, &
332	& Cd , Ch , Ce )
333	ELSE
334	!! If air temp. and spec. hum. are at a different height to wind :
335	CALL TURB_CORE_2Z(rn_zqt, rn_zu , zst , sf(jp_tair)%fnow, &
336	& zqsatw, sf(jp_humi)%fnow, wndm, &
337	& Cd , Ch , Ce , &
338	& zt_zu , zq_zu )
339	ENDIF
340	ELSE
341	!! If air temp. and spec. hum. are given at same height than wind (10m) :
342	!gm bug? at the compiling phase, add a copy in temporary arrays... ==> check perf
343	! CALL TURB_CORE_1Z( 10., zst (:,:), sf(jp_tair)%fnow(:,:), &
344	! & zqsatw(:,:), sf(jp_humi)%fnow(:,:), wndm(:,:), &
345	! & Cd (:,:), Ch (:,:), Ce (:,:) )
346	!gm bug
347	! ARPDBG - this won't compile with gfortran. Fix but check performance
348	! as per comment above.
349	CALL TURB_CORE_1Z( 10., zst , sf(jp_tair)%fnow(:,:,1), &
350	& zqsatw, sf(jp_humi)%fnow(:,:,1), wndm, &
351	& Cd , Ch , Ce )
352	ENDIF
353
354	! ... tau module, i and j component
355	DO jj = 1, jpj
356	DO ji = 1, jpi
357	zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj)
358	taum (ji,jj) = zztmp * wndm (ji,jj)
359	zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
360	zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
361	END DO
362	END DO
363
364	! ... add the HF tau contribution to the wind stress module?
365	IF( lhftau ) THEN
366	!CDIR COLLAPSE
367	taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
368	ENDIF
369	CALL iom_put( "taum_oce", taum ) ! output wind stress module
370
371	! ... utau, vtau at U- and V_points, resp.
372	! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
373	DO jj = 1, jpjm1
374	DO ji = 1, fs_jpim1
375	utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) )
376	vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) )
377	END DO
378	END DO
379	CALL lbc_lnk( utau(:,:), 'U', -1. )
380	CALL lbc_lnk( vtau(:,:), 'V', -1. )
381
382	! Turbulent fluxes over ocean
383	! -----------------------------
384	IF( ln_2m .OR. ( ln_bulk2z .AND. rn_zqt /= rn_zu ) ) THEN
385	! Values of temp. and hum. adjusted to height of wind must be used
386	zevap(:,:) = rn_efac * MAX( 0.e0, rhoa Ce(:,:)( zqsatw(:,:) - zq_zu(:,:) ) * wndm(:,:) ) ! Evaporation
387	zqsb (:,:) = rhoacpaCh(:,:)( zst (:,:) - zt_zu(:,:) ) wndm(:,:) ! Sensible Heat
388	ELSE
389	!CDIR COLLAPSE
390	zevap(:,:) = rn_efac * MAX( 0.e0, rhoa Ce(:,:)( zqsatw(:,:) - sf(jp_humi)%fnow(:,:,1) ) * wndm(:,:) ) ! Evaporation
391	!CDIR COLLAPSE
392	zqsb (:,:) = rhoacpaCh(:,:)( zst (:,:) - sf(jp_tair)%fnow(:,:,1) ) wndm(:,:) ! Sensible Heat
393	ENDIF
394	!CDIR COLLAPSE
395	zqla (:,:) = Lv * zevap(:,:) ! Latent Heat
396
397	IF(ln_ctl) THEN
398	CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_core: zqla : ', tab2d_2=Ce , clinfo2=' Ce : ' )
399	CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce_core: zqsb : ', tab2d_2=Ch , clinfo2=' Ch : ' )
400	CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_core: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' )
401	CALL prt_ctl( tab2d_1=zqsatw, clinfo1=' blk_oce_core: zqsatw : ', tab2d_2=zst, clinfo2=' zst : ' )
402	CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_core: utau : ', mask1=umask, &
403	& tab2d_2=vtau , clinfo2= ' vtau : ' , mask2=vmask )
404	CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_core: wndm : ')
405	CALL prt_ctl( tab2d_1=zst , clinfo1=' blk_oce_core: zst : ')
406	ENDIF
407
408	! ----------------------------------------------------------------------------- !
409	! III Total FLUXES !
410	! ----------------------------------------------------------------------------- !
411
412	!CDIR COLLAPSE
413	emp (:,:) = ( zevap(:,:) & ! mass flux (evap. - precip.)
414	& - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1)
415	!CDIR COLLAPSE
416	qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar flux
417	& - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus & ! remove latent melting heat for solid precip
418	& - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST
419	& + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac & ! add liquid precip heat content at Tair
420	& * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp &
421	& + sf(jp_snow)%fnow(:,:,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)
422	& * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic
423	!
424	CALL iom_put( "qlw_oce", zqlw ) ! output downward longwave heat over the ocean
425	CALL iom_put( "qsb_oce", - zqsb ) ! output downward sensible heat over the ocean
426	CALL iom_put( "qla_oce", - zqla ) ! output downward latent heat over the ocean
427	CALL iom_put( "qhc_oce", qns-zqlw+zqsb+zqla ) ! output downward heat content of E-P over the ocean
428	CALL iom_put( "qns_oce", qns ) ! output downward non solar heat over the ocean
429	!
430	IF(ln_ctl) THEN
431	CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce_core: zqsb : ', tab2d_2=zqlw , clinfo2=' zqlw : ')
432	CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_core: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ')
433	CALL prt_ctl(tab2d_1=pst , clinfo1=' blk_oce_core: pst : ', tab2d_2=emp , clinfo2=' emp : ')
434	CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce_core: utau : ', mask1=umask, &
435	& tab2d_2=vtau , clinfo2= ' vtau : ' , mask2=vmask )
436	ENDIF
437	!
438	CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap )
439	CALL wrk_dealloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )
440	!
441	IF( nn_timing == 1 ) CALL timing_stop('blk_oce_core')
442	!
443	END SUBROUTINE blk_oce_core
444
445	SUBROUTINE blk_bio_meanqsr
446	!!---------------------------------------------------------------------
447	!! *** ROUTINE blk_bio_meanqsr
448	!!
449	!! ** Purpose : provide daily qsr_mean for PISCES when
450	!! analytic diurnal cycle is applied in physic
451	!!
452	!! ** Method : add part where there is no ice
453	!!
454	!!---------------------------------------------------------------------
455	IF( nn_timing == 1 ) CALL timing_start('blk_bio_meanqsr')
456
457	qsr_mean(:,:) = (1. - albo ) * sf(jp_qsr)%fnow(:,:,1)
458
459	IF( nn_timing == 1 ) CALL timing_stop('blk_bio_meanqsr')
460
461	END SUBROUTINE blk_bio_meanqsr
462
463
464	SUBROUTINE blk_ice_meanqsr(palb,p_qsr_mean,pdim)
465	!!---------------------------------------------------------------------
466	!!
467	!! ** Purpose : provide the daily qsr_mean over sea_ice for PISCES when
468	!! analytic diurnal cycle is applied in physic
469	!!
470	!! ** Method : compute qsr
471	!!
472	!!---------------------------------------------------------------------
473	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb ! ice albedo (clear sky) (alb_ice_cs) [%]
474	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_qsr_mean ! solar heat flux over ice (T-point) [W/m2]
475	INTEGER , INTENT(in ) :: pdim ! number of ice categories
476	!!
477	INTEGER :: ijpl ! number of ice categories (size of 3rd dim of input arrays)
478	INTEGER :: ji, jj, jl ! dummy loop indices
479	REAL(wp) :: zztmp ! temporary variable
480	!!---------------------------------------------------------------------
481	IF( nn_timing == 1 ) CALL timing_start('blk_ice_meanqsr')
482	!
483	ijpl = pdim ! number of ice categories
484	zztmp = 1. / ( 1. - albo )
485	! ! ========================== !
486	DO jl = 1, ijpl ! Loop over ice categories !
487	! ! ========================== !
488	DO jj = 1 , jpj
489	DO ji = 1, jpi
490	p_qsr_mean(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr_mean(ji,jj)
491	END DO
492	END DO
493	END DO
494	!
495	IF( nn_timing == 1 ) CALL timing_stop('blk_ice_meanqsr')
496	!
497	END SUBROUTINE blk_ice_meanqsr
498
499
500	SUBROUTINE blk_ice_core( pst , pui , pvi , palb , &
501	& p_taui, p_tauj, p_qns , p_qsr, &
502	& p_qla , p_dqns, p_dqla, &
503	& p_tpr , p_spr , &
504	& p_fr1 , p_fr2 , cd_grid, pdim )
505	!!---------------------------------------------------------------------
506	!! * ROUTINE blk_ice_core *
507	!!
508	!! ** Purpose : provide the surface boundary condition over sea-ice
509	!!
510	!! ** Method : compute momentum, heat and freshwater exchanged
511	!! between atmosphere and sea-ice using CORE bulk
512	!! formulea, ice variables and read atmmospheric fields.
513	!! NB: ice drag coefficient is assumed to be a constant
514	!!
515	!! caution : the net upward water flux has with mm/day unit
516	!!---------------------------------------------------------------------
517	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pst ! ice surface temperature (>0, =rt0 over land) [Kelvin]
518	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pui ! ice surface velocity (i- and i- components [m/s]
519	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvi ! at I-point (B-grid) or U & V-point (C-grid)
520	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb ! ice albedo (clear sky) (alb_ice_cs) [%]
521	REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_taui ! i- & j-components of surface ice stress [N/m2]
522	REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_tauj ! at I-point (B-grid) or U & V-point (C-grid)
523	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_qns ! non solar heat flux over ice (T-point) [W/m2]
524	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_qsr ! solar heat flux over ice (T-point) [W/m2]
525	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_qla ! latent heat flux over ice (T-point) [W/m2]
526	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_dqns ! non solar heat sensistivity (T-point) [W/m2]
527	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_dqla ! latent heat sensistivity (T-point) [W/m2]
528	REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_tpr ! total precipitation (T-point) [Kg/m2/s]
529	REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_spr ! solid precipitation (T-point) [Kg/m2/s]
530	REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_fr1 ! 1sr fraction of qsr penetration in ice (T-point) [%]
531	REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_fr2 ! 2nd fraction of qsr penetration in ice (T-point) [%]
532	CHARACTER(len=1) , INTENT(in ) :: cd_grid ! ice grid ( C or B-grid)
533	INTEGER , INTENT(in ) :: pdim ! number of ice categories
534	!!
535	INTEGER :: ji, jj, jl ! dummy loop indices
536	INTEGER :: ijpl ! number of ice categories (size of 3rd dim of input arrays)
537	REAL(wp) :: zst2, zst3
538	REAL(wp) :: zcoef_wnorm, zcoef_wnorm2, zcoef_dqlw, zcoef_dqla, zcoef_dqsb
539	REAL(wp) :: zztmp ! temporary variable
540	REAL(wp) :: zcoef_frca ! fractional cloud amount
541	REAL(wp) :: zwnorm_f, zwndi_f , zwndj_f ! relative wind module and components at F-point
542	REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point
543	!!
544	REAL(wp), DIMENSION(:,:) , POINTER :: z_wnds_t ! wind speed ( = \| U10m - U_ice \| ) at T-point
545	REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice
546	REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice
547	REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqlw ! long wave heat sensitivity over ice
548	REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqsb ! sensible heat sensitivity over ice
549	!!---------------------------------------------------------------------
550	!
551	IF( nn_timing == 1 ) CALL timing_start('blk_ice_core')
552	!
553	CALL wrk_alloc( jpi,jpj, z_wnds_t )
554	CALL wrk_alloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb )
555
556	ijpl = pdim ! number of ice categories
557
558	! local scalars ( place there for vector optimisation purposes)
559	zcoef_wnorm = rhoa * Cice
560	zcoef_wnorm2 = rhoa * Cice * 0.5
561	zcoef_dqlw = 4.0 * 0.95 * Stef
562	zcoef_dqla = -Ls * Cice * 11637800. * (-5897.8)
563	zcoef_dqsb = rhoa * cpa * Cice
564	zcoef_frca = 1.0 - 0.3
565
566	!!gm brutal....
567	z_wnds_t(:,:) = 0.e0
568	p_taui (:,:) = 0.e0
569	p_tauj (:,:) = 0.e0
570	!!gm end
571
572	#if defined key_lim3
573	tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) ! LIM3: make Tair available in sea-ice. WARNING allocated after call to ice_init
574	#endif
575	! ----------------------------------------------------------------------------- !
576	! Wind components and module relative to the moving ocean ( U10m - U_ice ) !
577	! ----------------------------------------------------------------------------- !
578	SELECT CASE( cd_grid )
579	CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation)
580	! and scalar wind at T-point ( = \| U10m - U_ice \| ) (masked)
581	!CDIR NOVERRCHK
582	DO jj = 2, jpjm1
583	DO ji = 2, jpim1 ! B grid : NO vector opt
584	! ... scalar wind at I-point (fld being at T-point)
585	zwndi_f = 0.25 * ( sf(jp_wndi)%fnow(ji-1,jj ,1) + sf(jp_wndi)%fnow(ji ,jj ,1) &
586	& + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji ,jj-1,1) ) - rn_vfac * pui(ji,jj)
587	zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) &
588	& + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * pvi(ji,jj)
589	zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
590	! ... ice stress at I-point
591	p_taui(ji,jj) = zwnorm_f * zwndi_f
592	p_tauj(ji,jj) = zwnorm_f * zwndj_f
593	! ... scalar wind at T-point (fld being at T-point)
594	zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( pui(ji,jj+1) + pui(ji+1,jj+1) &
595	& + pui(ji,jj ) + pui(ji+1,jj ) )
596	zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( pvi(ji,jj+1) + pvi(ji+1,jj+1) &
597	& + pvi(ji,jj ) + pvi(ji+1,jj ) )
598	z_wnds_t(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
599	END DO
600	END DO
601	CALL lbc_lnk( p_taui , 'I', -1. )
602	CALL lbc_lnk( p_tauj , 'I', -1. )
603	CALL lbc_lnk( z_wnds_t, 'T', 1. )
604	!
605	CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean)
606	#if defined key_vectopt_loop
607	!CDIR COLLAPSE
608	#endif
609	DO jj = 2, jpj
610	DO ji = fs_2, jpi ! vect. opt.
611	zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pui(ji-1,jj ) + pui(ji,jj) ) )
612	zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pvi(ji ,jj-1) + pvi(ji,jj) ) )
613	z_wnds_t(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
614	END DO
615	END DO
616	#if defined key_vectopt_loop
617	!CDIR COLLAPSE
618	#endif
619	DO jj = 2, jpjm1
620	DO ji = fs_2, fs_jpim1 ! vect. opt.
621	p_taui(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji+1,jj ) + z_wnds_t(ji,jj) ) &
622	& * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * pui(ji,jj) )
623	p_tauj(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji,jj+1 ) + z_wnds_t(ji,jj) ) &
624	& * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * pvi(ji,jj) )
625	END DO
626	END DO
627	CALL lbc_lnk( p_taui , 'U', -1. )
628	CALL lbc_lnk( p_tauj , 'V', -1. )
629	CALL lbc_lnk( z_wnds_t, 'T', 1. )
630	!
631	END SELECT
632
633	zztmp = 1. / ( 1. - albo )
634	! ! ========================== !
635	DO jl = 1, ijpl ! Loop over ice categories !
636	! ! ========================== !
637	!CDIR NOVERRCHK
638	!CDIR COLLAPSE
639	DO jj = 1 , jpj
640	!CDIR NOVERRCHK
641	DO ji = 1, jpi
642	! ----------------------------!
643	! I Radiative FLUXES !
644	! ----------------------------!
645	zst2 = pst(ji,jj,jl) * pst(ji,jj,jl)
646	zst3 = pst(ji,jj,jl) * zst2
647	! Short Wave (sw)
648	p_qsr(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
649	! Long Wave (lw)
650	! iovino
651	IF( ff(ji,jj) .GT. 0._wp ) THEN
652	z_qlw(ji,jj,jl) = ( 0.95 * sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
653	ELSE
654	z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
655	ENDIF
656	! lw sensitivity
657	z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
658
659	! ----------------------------!
660	! II Turbulent FLUXES !
661	! ----------------------------!
662
663	! ... turbulent heat fluxes
664	! Sensible Heat
665	z_qsb(ji,jj,jl) = rhoa * cpa * Cice * z_wnds_t(ji,jj) * ( pst(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
666	! Latent Heat
667	p_qla(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * Cice * z_wnds_t(ji,jj) &
668	& * ( 11637800. * EXP( -5897.8 / pst(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1) ) )
669	! Latent heat sensitivity for ice (Dqla/Dt)
670	p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) )
671	! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
672	z_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj)
673
674	! ----------------------------!
675	! III Total FLUXES !
676	! ----------------------------!
677	! Downward Non Solar flux
678	p_qns (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - p_qla (ji,jj,jl)
679	! Total non solar heat flux sensitivity for ice
680	p_dqns(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + p_dqla(ji,jj,jl) )
681	END DO
682	!
683	END DO
684	!
685	END DO
686	!
687	!--------------------------------------------------------------------
688	! FRACTIONs of net shortwave radiation which is not absorbed in the
689	! thin surface layer and penetrates inside the ice cover
690	! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
691
692	!CDIR COLLAPSE
693	p_fr1(:,:) = ( 0.18 * ( 1.0 - zcoef_frca ) + 0.35 * zcoef_frca )
694	!CDIR COLLAPSE
695	p_fr2(:,:) = ( 0.82 * ( 1.0 - zcoef_frca ) + 0.65 * zcoef_frca )
696
697	!CDIR COLLAPSE
698	p_tpr(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac ! total precipitation [kg/m2/s]
699	!CDIR COLLAPSE
700	p_spr(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac ! solid precipitation [kg/m2/s]
701	CALL iom_put( 'snowpre', p_spr * 86400. ) ! Snow precipitation
702	CALL iom_put( 'precip', p_tpr * 86400. ) ! Total precipitation
703	!
704	IF(ln_ctl) THEN
705	CALL prt_ctl(tab3d_1=p_qla , clinfo1=' blk_ice_core: p_qla : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=ijpl)
706	CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice_core: z_qlw : ', tab3d_2=p_dqla , clinfo2=' p_dqla : ', kdim=ijpl)
707	CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice_core: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=ijpl)
708	CALL prt_ctl(tab3d_1=p_dqns , clinfo1=' blk_ice_core: p_dqns : ', tab3d_2=p_qsr , clinfo2=' p_qsr : ', kdim=ijpl)
709	CALL prt_ctl(tab3d_1=pst , clinfo1=' blk_ice_core: pst : ', tab3d_2=p_qns , clinfo2=' p_qns : ', kdim=ijpl)
710	CALL prt_ctl(tab2d_1=p_tpr , clinfo1=' blk_ice_core: p_tpr : ', tab2d_2=p_spr , clinfo2=' p_spr : ')
711	CALL prt_ctl(tab2d_1=p_taui , clinfo1=' blk_ice_core: p_taui : ', tab2d_2=p_tauj , clinfo2=' p_tauj : ')
712	CALL prt_ctl(tab2d_1=z_wnds_t, clinfo1=' blk_ice_core: z_wnds_t : ')
713	ENDIF
714
715	CALL wrk_dealloc( jpi,jpj, z_wnds_t )
716	CALL wrk_dealloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb )
717	!
718	IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core')
719	!
720	END SUBROUTINE blk_ice_core
721
722
723	SUBROUTINE TURB_CORE_1Z(zu, sst, T_a, q_sat, q_a, &
724	& dU , Cd , Ch , Ce )
725	!!----------------------------------------------------------------------
726	!! * ROUTINE turb_core *
727	!!
728	!! ** Purpose : Computes turbulent transfert coefficients of surface
729	!! fluxes according to Large & Yeager (2004)
730	!!
731	!! ** Method : I N E R T I A L D I S S I P A T I O N M E T H O D
732	!! Momentum, Latent and sensible heat exchange coefficients
733	!! Caution: this procedure should only be used in cases when air
734	!! temperature (T_air), air specific humidity (q_air) and wind (dU)
735	!! are provided at the same height 'zzu'!
736	!!
737	!! References : Large & Yeager, 2004 : ???
738	!!----------------------------------------------------------------------
739	REAL(wp) , INTENT(in ) :: zu ! altitude of wind measurement [m]
740	REAL(wp), DIMENSION(:,:), INTENT(in ) :: sst ! sea surface temperature [Kelvin]
741	REAL(wp), DIMENSION(:,:), INTENT(in ) :: T_a ! potential air temperature [Kelvin]
742	REAL(wp), DIMENSION(:,:), INTENT(in ) :: q_sat ! sea surface specific humidity [kg/kg]
743	REAL(wp), DIMENSION(:,:), INTENT(in ) :: q_a ! specific air humidity [kg/kg]
744	REAL(wp), DIMENSION(:,:), INTENT(in ) :: dU ! wind module \|U(zu)-U(0)\| [m/s]
745	REAL(wp), DIMENSION(:,:), INTENT( out) :: Cd ! transfert coefficient for momentum (tau)
746	REAL(wp), DIMENSION(:,:), INTENT( out) :: Ch ! transfert coefficient for temperature (Q_sens)
747	REAL(wp), DIMENSION(:,:), INTENT( out) :: Ce ! transfert coefficient for evaporation (Q_lat)
748	!!
749	INTEGER :: j_itt
750	INTEGER , PARAMETER :: nb_itt = 3
751	REAL(wp), PARAMETER :: grav = 9.8 ! gravity
752	REAL(wp), PARAMETER :: kappa = 0.4 ! von Karman s constant
753
754	REAL(wp), DIMENSION(:,:), POINTER :: dU10 ! dU [m/s]
755	REAL(wp), DIMENSION(:,:), POINTER :: dT ! air/sea temperature differeence [K]
756	REAL(wp), DIMENSION(:,:), POINTER :: dq ! air/sea humidity difference [K]
757	REAL(wp), DIMENSION(:,:), POINTER :: Cd_n10 ! 10m neutral drag coefficient
758	REAL(wp), DIMENSION(:,:), POINTER :: Ce_n10 ! 10m neutral latent coefficient
759	REAL(wp), DIMENSION(:,:), POINTER :: Ch_n10 ! 10m neutral sensible coefficient
760	REAL(wp), DIMENSION(:,:), POINTER :: sqrt_Cd_n10 ! root square of Cd_n10
761	REAL(wp), DIMENSION(:,:), POINTER :: sqrt_Cd ! root square of Cd
762	REAL(wp), DIMENSION(:,:), POINTER :: T_vpot ! virtual potential temperature [K]
763	REAL(wp), DIMENSION(:,:), POINTER :: T_star ! turbulent scale of tem. fluct.
764	REAL(wp), DIMENSION(:,:), POINTER :: q_star ! turbulent humidity of temp. fluct.
765	REAL(wp), DIMENSION(:,:), POINTER :: U_star ! turb. scale of velocity fluct.
766	REAL(wp), DIMENSION(:,:), POINTER :: L ! Monin-Obukov length [m]
767	REAL(wp), DIMENSION(:,:), POINTER :: zeta ! stability parameter at height zu
768	REAL(wp), DIMENSION(:,:), POINTER :: U_n10 ! neutral wind velocity at 10m [m]
769	REAL(wp), DIMENSION(:,:), POINTER :: xlogt, xct, zpsi_h, zpsi_m
770
771	INTEGER , DIMENSION(:,:), POINTER :: stab ! 1st guess stability test integer
772	!!----------------------------------------------------------------------
773	!
774	IF( nn_timing == 1 ) CALL timing_start('TURB_CORE_1Z')
775	!
776	CALL wrk_alloc( jpi,jpj, stab ) ! integer
777	CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
778	CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m )
779
780	!! * Start
781	!! Air/sea differences
782	dU10 = max(0.5, dU) ! we don't want to fall under 0.5 m/s
783	dT = T_a - sst ! assuming that T_a is allready the potential temp. at zzu
784	dq = q_a - q_sat
785	!!
786	!! Virtual potential temperature
787	T_vpot = T_a(1. + 0.608q_a)
788	!!
789	!! Neutral Drag Coefficient
790	stab = 0.5 + sign(0.5,dT) ! stable : stab = 1 ; unstable : stab = 0
791	IF ( ln_cdgw ) THEN
792	cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
793	Cd_n10(:,:) = cdn_wave
794	ELSE
795	Cd_n10 = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a)
796	ENDIF
797	sqrt_Cd_n10 = sqrt(Cd_n10)
798	Ce_n10 = 1.e-3 * ( 34.6 * sqrt_Cd_n10 ) ! L & Y eq. (6b)
799	Ch_n10 = 1.e-3sqrt_Cd_n10(18.stab + 32.7(1.-stab)) ! L & Y eq. (6c), (6d)
800	!!
801	!! Initializing transfert coefficients with their first guess neutral equivalents :
802	Cd = Cd_n10 ; Ce = Ce_n10 ; Ch = Ch_n10 ; sqrt_Cd = sqrt(Cd)
803
804	!! * Now starting iteration loop
805	DO j_itt=1, nb_itt
806	!! Turbulent scales :
807	U_star = sqrt_Cd*dU10 ! L & Y eq. (7a)
808	T_star = Ch/sqrt_Cd*dT ! L & Y eq. (7b)
809	q_star = Ce/sqrt_Cd*dq ! L & Y eq. (7c)
810
811	!! Estimate the Monin-Obukov length :
812	L = (U_star*2)/( kappagrav*(T_star/T_vpot + q_star/(q_a + 1./0.608)) )
813
814	!! Stability parameters :
815	zeta = zu/L ; zeta = sign( min(abs(zeta),10.0), zeta )
816	zpsi_h = psi_h(zeta)
817	zpsi_m = psi_m(zeta)
818
819	IF ( ln_cdgw ) THEN
820	sqrt_Cd=kappa/((kappa/sqrt_Cd_n10) - zpsi_m) ; Cd=sqrt_Cd*sqrt_Cd;
821	ELSE
822	!! Shifting the wind speed to 10m and neutral stability :
823	U_n10 = dU101./(1. + sqrt_Cd_n10/kappa(log(zu/10.) - zpsi_m)) ! L & Y eq. (9a)
824
825	!! Updating the neutral 10m transfer coefficients :
826	Cd_n10 = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a)
827	sqrt_Cd_n10 = sqrt(Cd_n10)
828	Ce_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b)
829	stab = 0.5 + sign(0.5,zeta)
830	Ch_n10 = 1.e-3sqrt_Cd_n10(18.stab + 32.7(1.-stab)) ! L & Y eq. (6c), (6d)
831
832	!! Shifting the neutral 10m transfer coefficients to ( zu , zeta ) :
833	!!
834	xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10) - zpsi_m)
835	Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd)
836	ENDIF
837	!!
838	xlogt = log(zu/10.) - zpsi_h
839	!!
840	xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10
841	Ch = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct
842	!!
843	xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10
844	Ce = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct
845	!!
846	END DO
847	!!
848	CALL wrk_dealloc( jpi,jpj, stab ) ! integer
849	CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
850	CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m )
851	!
852	IF( nn_timing == 1 ) CALL timing_stop('TURB_CORE_1Z')
853	!
854	END SUBROUTINE TURB_CORE_1Z
855
856
857	SUBROUTINE TURB_CORE_2Z(zt, zu, sst, T_zt, q_sat, q_zt, dU, Cd, Ch, Ce, T_zu, q_zu)
858	!!----------------------------------------------------------------------
859	!! * ROUTINE turb_core *
860	!!
861	!! ** Purpose : Computes turbulent transfert coefficients of surface
862	!! fluxes according to Large & Yeager (2004).
863	!!
864	!! ** Method : I N E R T I A L D I S S I P A T I O N M E T H O D
865	!! Momentum, Latent and sensible heat exchange coefficients
866	!! Caution: this procedure should only be used in cases when air
867	!! temperature (T_air) and air specific humidity (q_air) are at a
868	!! different height to wind (dU).
869	!!
870	!! References : Large & Yeager, 2004 : ???
871	!!----------------------------------------------------------------------
872	REAL(wp), INTENT(in ) :: zt ! height for T_zt and q_zt [m]
873	REAL(wp), INTENT(in ) :: zu ! height for dU [m]
874	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst ! sea surface temperature [Kelvin]
875	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: T_zt ! potential air temperature [Kelvin]
876	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_sat ! sea surface specific humidity [kg/kg]
877	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity [kg/kg]
878	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: dU ! relative wind module \|U(zu)-U(0)\| [m/s]
879	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau)
880	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens)
881	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ce ! transfert coefficient for evaporation (Q_lat)
882	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: T_zu ! air temp. shifted at zu [K]
883	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. hum. shifted at zu [kg/kg]
884
885	INTEGER :: j_itt
886	INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations
887	REAL(wp), PARAMETER :: grav = 9.8 ! gravity
888	REAL(wp), PARAMETER :: kappa = 0.4 ! von Karman's constant
889
890	REAL(wp), DIMENSION(:,:), POINTER :: dU10 ! dU [m/s]
891	REAL(wp), DIMENSION(:,:), POINTER :: dT ! air/sea temperature differeence [K]
892	REAL(wp), DIMENSION(:,:), POINTER :: dq ! air/sea humidity difference [K]
893	REAL(wp), DIMENSION(:,:), POINTER :: Cd_n10 ! 10m neutral drag coefficient
894	REAL(wp), DIMENSION(:,:), POINTER :: Ce_n10 ! 10m neutral latent coefficient
895	REAL(wp), DIMENSION(:,:), POINTER :: Ch_n10 ! 10m neutral sensible coefficient
896	REAL(wp), DIMENSION(:,:), POINTER :: sqrt_Cd_n10 ! root square of Cd_n10
897	REAL(wp), DIMENSION(:,:), POINTER :: sqrt_Cd ! root square of Cd
898	REAL(wp), DIMENSION(:,:), POINTER :: T_vpot ! virtual potential temperature [K]
899	REAL(wp), DIMENSION(:,:), POINTER :: T_star ! turbulent scale of tem. fluct.
900	REAL(wp), DIMENSION(:,:), POINTER :: q_star ! turbulent humidity of temp. fluct.
901	REAL(wp), DIMENSION(:,:), POINTER :: U_star ! turb. scale of velocity fluct.
902	REAL(wp), DIMENSION(:,:), POINTER :: L ! Monin-Obukov length [m]
903	REAL(wp), DIMENSION(:,:), POINTER :: zeta_u ! stability parameter at height zu
904	REAL(wp), DIMENSION(:,:), POINTER :: zeta_t ! stability parameter at height zt
905	REAL(wp), DIMENSION(:,:), POINTER :: U_n10 ! neutral wind velocity at 10m [m]
906	REAL(wp), DIMENSION(:,:), POINTER :: xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m
907
908	INTEGER , DIMENSION(:,:), POINTER :: stab ! 1st stability test integer
909	!!----------------------------------------------------------------------
910	!
911	IF( nn_timing == 1 ) CALL timing_start('TURB_CORE_2Z')
912	!
913	CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
914	CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 )
915	CALL wrk_alloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m )
916	CALL wrk_alloc( jpi,jpj, stab ) ! interger
917
918	!! Initial air/sea differences
919	dU10 = max(0.5, dU) ! we don't want to fall under 0.5 m/s
920	dT = T_zt - sst
921	dq = q_zt - q_sat
922
923	!! Neutral Drag Coefficient :
924	stab = 0.5 + sign(0.5,dT) ! stab = 1 if dT > 0 -> STABLE
925	IF( ln_cdgw ) THEN
926	cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
927	Cd_n10(:,:) = cdn_wave
928	ELSE
929	Cd_n10 = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 )
930	ENDIF
931	sqrt_Cd_n10 = sqrt(Cd_n10)
932	Ce_n10 = 1.e-3( 34.6 sqrt_Cd_n10 )
933	Ch_n10 = 1.e-3sqrt_Cd_n10(18.stab + 32.7(1. - stab))
934
935	!! Initializing transf. coeff. with their first guess neutral equivalents :
936	Cd = Cd_n10 ; Ce = Ce_n10 ; Ch = Ch_n10 ; sqrt_Cd = sqrt(Cd)
937
938	!! Initializing z_u values with z_t values :
939	T_zu = T_zt ; q_zu = q_zt
940
941	!! * Now starting iteration loop
942	DO j_itt=1, nb_itt
943	dT = T_zu - sst ; dq = q_zu - q_sat ! Updating air/sea differences
944	T_vpot = T_zu(1. + 0.608q_zu) ! Updating virtual potential temperature at zu
945	U_star = sqrt_Cd*dU10 ! Updating turbulent scales : (L & Y eq. (7))
946	T_star = Ch/sqrt_Cd*dT !
947	q_star = Ce/sqrt_Cd*dq !
948	!!
949	L = (U_star*U_star) & ! Estimate the Monin-Obukov length at height zu
950	& / (kappagrav/T_vpot(T_star(1.+0.608q_zu) + 0.608T_zuq_star))
951	!! Stability parameters :
952	zeta_u = zu/L ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u )
953	zeta_t = zt/L ; zeta_t = sign( min(abs(zeta_t),10.0), zeta_t )
954	zpsi_hu = psi_h(zeta_u)
955	zpsi_ht = psi_h(zeta_t)
956	zpsi_m = psi_m(zeta_u)
957	!!
958	!! Shifting the wind speed to 10m and neutral stability : (L & Y eq.(9a))
959	! U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u)))
960	U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m))
961	!!
962	!! Shifting temperature and humidity at zu : (L & Y eq. (9b-9c))
963	! T_zu = T_zt - T_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t))
964	T_zu = T_zt - T_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht)
965	! q_zu = q_zt - q_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t))
966	q_zu = q_zt - q_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht)
967	!!
968	!! q_zu cannot have a negative value : forcing 0
969	stab = 0.5 + sign(0.5,q_zu) ; q_zu = stab*q_zu
970	!!
971	IF( ln_cdgw ) THEN
972	sqrt_Cd=kappa/((kappa/sqrt_Cd_n10) - zpsi_m) ; Cd=sqrt_Cd*sqrt_Cd;
973	ELSE
974	!! Updating the neutral 10m transfer coefficients :
975	Cd_n10 = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a)
976	sqrt_Cd_n10 = sqrt(Cd_n10)
977	Ce_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b)
978	stab = 0.5 + sign(0.5,zeta_u)
979	Ch_n10 = 1.e-3sqrt_Cd_n10(18.stab + 32.7(1.-stab)) ! L & Y eq. (6c-6d)
980	!!
981	!!
982	!! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) :
983	xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m) ! L & Y eq. (10a)
984	Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd)
985	ENDIF
986	!!
987	xlogt = log(zu/10.) - zpsi_hu
988	!!
989	xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10 ! L & Y eq. (10b)
990	Ch = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct
991	!!
992	xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10 ! L & Y eq. (10c)
993	Ce = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct
994	!!
995	!!
996	END DO
997	!!
998	CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
999	CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 )
1000	CALL wrk_dealloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m )
1001	CALL wrk_dealloc( jpi,jpj, stab ) ! interger
1002	!
1003	IF( nn_timing == 1 ) CALL timing_stop('TURB_CORE_2Z')
1004	!
1005	END SUBROUTINE TURB_CORE_2Z
1006
1007
1008	FUNCTION psi_m(zta) !! Psis, L & Y eq. (8c), (8d), (8e)
1009	!-------------------------------------------------------------------------------
1010	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta
1011
1012	REAL(wp), PARAMETER :: pi = 3.141592653589793_wp
1013	REAL(wp), DIMENSION(jpi,jpj) :: psi_m
1014	REAL(wp), DIMENSION(:,:), POINTER :: X2, X, stabit
1015	!-------------------------------------------------------------------------------
1016
1017	CALL wrk_alloc( jpi,jpj, X2, X, stabit )
1018
1019	X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.0) ; X = sqrt(X2)
1020	stabit = 0.5 + sign(0.5,zta)
1021	psi_m = -5.ztastabit & ! Stable
1022	& + (1. - stabit)(2log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable
1023
1024	CALL wrk_dealloc( jpi,jpj, X2, X, stabit )
1025	!
1026	END FUNCTION psi_m
1027
1028
1029	FUNCTION psi_h( zta ) !! Psis, L & Y eq. (8c), (8d), (8e)
1030	!-------------------------------------------------------------------------------
1031	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta
1032	!
1033	REAL(wp), DIMENSION(jpi,jpj) :: psi_h
1034	REAL(wp), DIMENSION(:,:), POINTER :: X2, X, stabit
1035	!-------------------------------------------------------------------------------
1036
1037	CALL wrk_alloc( jpi,jpj, X2, X, stabit )
1038
1039	X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.) ; X = sqrt(X2)
1040	stabit = 0.5 + sign(0.5,zta)
1041	psi_h = -5.ztastabit & ! Stable
1042	& + (1. - stabit)(2.log( (1. + X2)/2. )) ! Unstable
1043
1044	CALL wrk_dealloc( jpi,jpj, X2, X, stabit )
1045	!
1046	END FUNCTION psi_h
1047
1048	!!======================================================================
1049	END MODULE sbcblk_core

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

Download in other formats: