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_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 @ 4040

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

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

Download in other formats: