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

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

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

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

Download in other formats: