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.

flx_core.h90 in trunk/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/SBC/flx_core.h90 @ 623

Last change on this file since 623 was 606, checked in by opalod, 17 years ago

nemo_v2_update_002 : CT : when using CORE forcing:

add a the possibility to use 2m air temperature and specific humidity with 10m wind speed via ln_2m logical
add possibility to use Katabatic winds enhancement via ln_kata logical

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 53.9 KB

Line
1	!!----------------------------------------------------------------------
2	!! * flx_core.h90 *
3	!!----------------------------------------------------------------------
4	!! flx : update surface thermohaline fluxes from the NCAR data set
5	!! read in a NetCDF file
6	!!----------------------------------------------------------------------
7	!! * Modules used C A U T I O N already defined in flxmod.F90
8
9	!! * Shared module variables
10	REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: &
11	gsss, & !: SSS mean on nfbulk ocean time step
12	gu , gv , & !: (un,vn)(jk=1) mean over nfbulk ocean time step
13	! ! these variables are used for output in diawri
14	qlw_oce , & !: long wave flx over ocean
15	qla_oce , & !: latent heat flx over ocean
16	qsb_oce , & !: sensible heat flx over ocean
17	qlw_ice , & !: long wave flx over ice
18	qsb_ice !: sensible heat flx over ice
19
20	!! * Module variables
21	INTEGER, PARAMETER :: jpfile = 8 ! maximum number of files to read
22	CHARACTER (LEN=32), DIMENSION (jpfile) :: clvarname
23	CHARACTER (LEN=50), DIMENSION (jpfile) :: clname
24	CHARACTER (LEN=32) :: clvarkatax , clvarkatay ! variable name for katabatic masks
25	CHARACTER (LEN=256) :: clnamekata ! filename for katabatic masks
26
27	INTEGER :: isnap
28	INTEGER, DIMENSION(jpfile) :: &
29	numflxall, & ! logical units for surface fluxes data
30	nrecflx , nrecflx2 ! first and second record to be read in flux file
31	REAL(wp), DIMENSION(jpfile) :: &
32	freqh ! Frequency for each forcing file (hours)
33	! ! a negative value of -12 corresponds to monthly
34	REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! used for modif wind stress in the first sea points
35	& rmskkatax, rmskkatay !: mask ocean to increase wind stress in first sea points
36	REAL(wp), DIMENSION(jpi,jpj,jpfile,2) :: &
37	flxdta !: set of NCAR 6hourly/daily/monthly fluxes
38	LOGICAL :: &
39	& ln_2m = .FALSE., & !: logical flag for height of air temp. and hum.
40	& ln_kata = .FALSE. !: logical flag for Katabatic winds enhancement.
41	!!----------------------------------------------------------------------
42	!! History :
43	!! 9.0 ! 04-08 (U. Schweckendiek) Original code
44	!! ! 05-04 (L. Brodeau, A.M. Treguier) severals modifications:
45	!! ! - new bulk routine for efficiency
46	!! ! - WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files
47	!! ! - file names and file characteristics in namelist
48	!! ! - Implement reading of 6-hourly fields
49	!! ! 06-02 (S. Masson, G. Madec) IOM read + NEMO "style"
50	!! ! 07-06 (P. Mathiot, J-M Molines) Katabatic winds enhancement
51	!! ! 12-06 (L. Brodeau) handle both 2m and 10m input fields
52	!!----------------------------------------------------------------------
53	!! OPA 9.0 , LOCEAN-IPSL (2006)
54	!! $Header$
55	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
56	!!----------------------------------------------------------------------
57
58	CONTAINS
59
60	SUBROUTINE flx( kt )
61	!!---------------------------------------------------------------------
62	!! * ROUTINE flx *
63	!!
64	!! ** Purpose : provide the thermohaline fluxes (heat and freshwater)
65	!! and momentum fluxes (tau)
66	!! to the ocean at each time step.
67	!!
68	!! ** Method : Read NCAR data in a NetCDF file
69	!! (file names and frequency of inputs specified in namelist)
70	!! precipitation: total (rain+snow)
71	!! precipitation: snow only
72	!! u10,v10 -> scalar wind at 10m in m/s - ON 'T' GRID POINTS!!!
73	!! solar radiation (short wave) in W/m2
74	!! thermal radiation (long wave) in W/m2
75	!! specific humidity in %
76	!! temperature at 10m in degrees K
77	!!
78	!! ** Action :
79	!! call flx_blk_albedo to compute ocean and ice albedo
80	!! Calculates forcing fluxes to input into ice model
81	! or to be used directly for the ocean in ocesbc.
82	!!
83	!! ** Outputs
84	!! COMMENTS TO BE ADDED -units to be verified!
85	!! taux : zonal wind stress on "u" points (N/m2)
86	!! tauy : meridional wind stress on "v" points (N/m2)
87	!! qsr_oce : Solar flux over the ocean (W/m2)
88	!! qnsr_oce: longwave flux over the ocean (W/m2)
89	!! qsr_ice : solar flux over the ice (W/m2)
90	!! qnsr_ice: longwave flux over the ice (W/m2)
91	!! qla_ice : latent flux over sea-ice (W/m2)
92	!! dqns_ice: total heat fluxes sensitivity over ice (dQ/dT)
93	!! dqla_ice: latent heat flux sensitivity over ice (dQla/dT)
94	!! fr1_i0
95	!! fr2_i0
96	!! tprecip
97	!! sprecip
98	!! evap
99	!!
100	!! caution : now, in the opa global model, the net upward water flux is
101	!! ------- with mm/day unit,i.e. Kg/m2/s.
102	!!---------------------------------------------------------------------
103	!! * modules used
104	USE iom
105	USE par_oce
106	USE flx_oce
107	USE blk_oce ! bulk variable
108	USE taumod
109	USE bulk
110	USE phycst
111	USE lbclnk
112	USE dtatem ! FOR Ce = F(SST(levitus)):
113
114	!! * arguments
115	INTEGER, INTENT( in ) :: kt ! ocean time step
116
117	!! * Local declarations
118	INTEGER , PARAMETER :: jpmaxtime = 366*4 ! maximum time steps in file will be for a 6 hourly field
119	! and a leap year (necessary variable for flinopen??)
120	INTEGER :: irectot, irecflx
121	INTEGER :: ihour ! current hour of the day at which we read the forcings
122	INTEGER :: &
123	imois, imois2, & ! temporary integers
124	i15 , iman , inum , & ! " "
125	ifpr , jfpr , & ! frequency of index for print in prihre
126	jj , ji ! Loop indices
127	REAL(wp) :: &
128	zadatrjb, & ! date (noninteger day) at which we read the forcings
129	zxy , & ! scalar for temporal interpolation
130	zcof ! scalar
131	REAL(wp), DIMENSION(jpi,jpj, jpfile) :: &
132	flxnow ! input flux values at current time step
133	REAL(wp), DIMENSION(jpi,jpj) :: &
134	dqlw_ice , & ! long wave heat flx sensitivity over ice
135	dqsb_ice , & ! sensible heat flx sensitivity over ice
136	alb_oce_os, & ! albedo of the ocean under overcast sky
137	alb_ice_os, & ! albedo of the ice under overcast sky
138	alb_ice_cs, & ! albedo of ice under clear sky
139	alb_oce_cs, & ! albedo of the ocean under clear sky
140	zsst, & ! nfbulk : mean SST
141	zsss, & ! nfbulk : mean tn_ice(:,:,1)
142	zut, & ! nfbulk : mean U at T-point
143	zvt, & ! nfbulk : mean V at T-point
144	dUnormt, & ! scalar wind (norm) on T points
145	tauxt, tauyt, & ! wind stress computed at T-point
146	qsatw, & ! specific humidity at zSST
147	qsat, & ! specific humidity at zsss
148	Ch, & ! coefficient for sensible heat flux
149	Ce, & ! coefficient for latent heat flux
150	Cd, & ! coefficient for wind stress
151	zt_zu, & !: air temperature at wind speed height
152	zq_zu !: air spec. hum. at wind speed height
153	REAL(wp), PARAMETER :: &
154	rhoa = 1.22, & ! air density
155	cp = 1000.5, & ! specific heat of air
156	Lv = 2.5e6, & ! latent heat of vaporization
157	Ls = 2.839e6, & ! latent heat of sublimation
158	Stef = 5.67e-8, & ! Stefan Boltzmann constant
159	Cice = 1.63e-3 ! transfert coefficient over ice
160	REAL(wp), DIMENSION(jpi,jpj) :: &
161	catm1 !
162
163	NAMELIST /namcore/ clname, clvarname, freqh, ln_2m, ln_kata
164	!!---------------------------------------------------------------------
165
166	!! =============================================== !!
167	!! PART A: READING FLUX FILES WHEN NECESSARY !!
168	!! =============================================== !!
169
170	!--- calculation for monthly data
171	i15 = INT( 2 * FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) )
172	iman = 12
173	imois = nmonth + i15 - 1
174	IF( imois == 0 ) imois = iman
175	imois2 = nmonth
176
177	!--- calculation for 6-hourly data
178	! we need the fraction of day, measured in hours at the date of forcings.
179	! This is the time step before the date we are calculating, zadatrj
180	! "adatrj" (real time in days (noninteger))
181	zadatrjb = adatrj - rdt / rday
182	ihour = INT( ( zadatrjb - INT(zadatrjb) ) * 24 )
183
184	! ! ------------------------ !
185	IF( kt == nit000 ) THEN ! first call kt=nit000 !
186	! ! ------------------------ !
187
188	!--- Initializes default values of file names, frequency of forcings
189	! and variable to be read in the files, before reading namelist
190	clname(1) = 'precip_core.nc' ; freqh(1) = -12 ; clvarname(1) = 'precip' ! monthly
191	clname(2) = 'u10_core.nc' ; freqh(2) = 24 ; clvarname(2) = 'u_10' ! 6 hourly
192	clname(3) = 'v10_core.nc' ; freqh(3) = 24 ; clvarname(3) = 'v_10' ! 6 hourly
193	clname(4) = 'q10_core.nc' ; freqh(4) = 24 ; clvarname(4) = 'q_10' ! 6 hourly
194	clname(5) = 'tot_solar_core.nc'; freqh(5) = 24 ; clvarname(5) = 'tot_solar' ! daily
195	clname(6) = 'therm_rad_core.nc'; freqh(6) = 24 ; clvarname(6) = 'therm_rad' ! daily
196	clname(7) = 'temp_10m_core.nc' ; freqh(7) = 24 ; clvarname(7) = 't_10' ! 6 hourly
197	clname(8) = 'snow_core.nc' ; freqh(8) = -12 ; clvarname(8) = 'snow' ! monthly
198
199	!--- Read Namelist namcore : OMIP/CORE
200	REWIND ( numnam )
201	READ ( numnam, namcore )
202	IF(lwp) THEN
203	WRITE(numout,*)' '
204	WRITE(numout,*)' flxmod/flx_core : global CORE fields in NetCDF format'
205	WRITE(numout,*)' ~~~~~~~~~~~~~~~ '
206	WRITE(numout,*) ' ln_2m = ', ln_2m
207	WRITE(numout,*) ' ln_kata = ', ln_kata
208	WRITE(numout,*) ' list of files and frequency (hour), or monthly (-12) '
209	DO ji = 1, jpfile
210	WRITE(numout,*) trim(clname(ji)), ' frequency:', freqh(ji)
211	END DO
212	ENDIF
213
214	!--- Initialization to zero
215	qsr_oce (:,:) = 0.e0
216	qsb_oce (:,:) = 0.e0
217	qla_oce (:,:) = 0.e0
218	qlw_oce (:,:) = 0.e0
219	qnsr_oce(:,:) = 0.e0
220	qsr_ice (:,:) = 0.e0
221	qnsr_ice(:,:) = 0.e0
222	qla_ice (:,:) = 0.e0
223	qlw_ice (:,:) = 0.e0
224	qsb_ice (:,:) = 0.e0
225
226	dqns_ice(:,:) = 0.e0
227	dqla_ice(:,:) = 0.e0
228	tprecip (:,:) = 0.e0
229	sprecip (:,:) = 0.e0
230	evap (:,:) = 0.e0
231
232	flxnow (:,:,:) = 1.e-15 !RB
233	flxdta (:,:,:,:) = 1.e-15 !RB
234
235	nrecflx (:) = 0 ! switch for reading flux data for each file
236	nrecflx2(:) = 0 ! switch for reading flux data for each file
237
238	!--- Open the files of the list
239	DO ji = 1, jpfile
240	CALL iom_open( clname(ji), numflxall(ji) )
241	END DO
242
243	ENDIF
244
245	! ! ------------------------ !
246	! ! Read files if required !
247	! ! ------------------------ !
248
249	!--- read data if necessary checks each file in turn
250
251	DO ji = 1, jpfile
252
253	IF ( freqh(ji) .GT.0 .AND. freqh(ji) .LT. 24 ) THEN !--- Case of 6-hourly flux data
254	!--- calculate current snapshot from hour of day
255	isnap = ihour / INT( freqh(ji) ) + 1
256	!--- reading is necessary when nrecflx(ji) differs from isnap
257	IF( nrecflx(ji) /= isnap ) THEN
258	nrecflx(ji) = isnap
259	irecflx = (nday_year-1)*24 / freqh(ji) + isnap
260	irectot = 365*24 / freqh(ji) ! this variable is not used by iom_get
261	IF(lwp) WRITE (numout,*)' CORE flux 6-h record :',irecflx, ' var:',trim(clvarname(ji))
262	CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), irecflx )
263	ENDIF
264
265	ELSE IF ( freqh(ji) .EQ. 24 ) THEN !--- Case of daily flux data
266
267	!--- reading is necessary when nrecflx(ji) differs from nday
268	IF( nrecflx(ji) /= nday ) THEN
269	nrecflx(ji) = nday !! remember present read day
270	irecflx = nday_year !!
271	irectot = 365 ! this variable is not used by iom_get
272	IF(lwp) WRITE (numout,*)'DAILY CORE flux record :',irecflx, ' var:',trim(clvarname(ji))
273	CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), irecflx )
274	ENDIF
275
276	ELSE IF ( freqh(ji) .EQ. -12 ) THEN !--- Case monthly data from CORE
277	! we read two months all the time although we
278	! could only read one and swap the arrays
279	IF( kt == nit000 .OR. imois /= nrecflx(ji) ) THEN
280	! calendar computation
281	! nrecflx number of the first file record used in the simulation
282	! nrecflx2 number of the last file record
283	nrecflx(ji) = imois
284	nrecflx2(ji) = nrecflx(ji)+1
285	nrecflx(ji) = MOD( nrecflx(ji), iman )
286	IF( nrecflx(ji) == 0 ) nrecflx(ji) = iman
287	nrecflx2(ji) = MOD( nrecflx2(ji), iman )
288	IF ( nrecflx2(ji) == 0 ) nrecflx2(ji) = iman
289	irectot = 12 ! this variable is not used by iom_get
290	IF(lwp) WRITE(numout,*) 'MONTHLY CORE flux record 1:', nrecflx(ji), &
291	& ' rec 2:', nrecflx2(ji), ' var:', trim(clvarname(ji))
292	CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), nrecflx (ji) )
293	CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,2), nrecflx2(ji) )
294	ENDIF
295	ENDIF
296	END DO ! end of loop over forcing files to verify if reading is necessary
297
298	!--- Printout if required
299	! IF( MOD( kt - 1 , nfbulk ) == 0 ) THEN
300	! printout at the initial time step only (unless you want to debug!!)
301	IF( kt == nit000 ) THEN
302	IF(lwp) THEN
303	WRITE(numout,*)
304	WRITE(numout,*) ' read daily and monthly CORE fluxes: ok'
305	WRITE(numout,*)
306	ifpr = INT(jpi/8) ; jfpr = INT(jpj/10)
307	DO ji = 1, jpfile
308	WRITE(numout,*) trim(clvarname(ji)),' day: ',ndastp
309	CALL prihre(flxdta(1,1,ji,1),jpi,jpj,1,jpi,ifpr,1,jpj,jfpr,0.,numout)
310	WRITE(numout,*)
311	END DO
312	ENDIF
313	ENDIF
314
315	! ! ------------------------ !
316	! ! Interpolates in time !
317	! ! ------------------------ !
318	! N.B. For now, only monthly values are interpolated, and they
319	! are interpolated to the current day, not to the time step.
320	!
321	DO ji = 1, jpfile
322	IF( freqh(ji) .EQ. -12 ) THEN
323	zxy = REAL(nday) / REAL( nobis(imois) ) + 0.5 - i15 !!! <== Caution nobis hard coded !!!
324	flxnow(:,:,ji) = ( ( 1. - zxy) * flxdta(:,:,ji,1) + zxy * flxdta(:,:,ji,2) )
325	ELSE
326	flxnow(:,:,ji) = flxdta(:,:,ji,1)
327	ENDIF
328	END DO
329	! JMM : add vatm needed in tracer routines ==> wind module ???
330	vatm(:,:) = SQRT( flxnow(:,:,2) * flxnow(:,:,2) + flxnow(:,:,3) * flxnow(:,:,3) )
331
332
333	!! =============================================== !!
334	!! PART B: CORE BULK CALCULATION !!
335	!! =============================================== !!
336
337	! for other forcing cases this is done in modules bulk.F90 and flxblk
338
339	! ! ------------------------ !
340	! ! Bulk initialisation !
341	! ! ------------------------ !
342
343	! code part from bulk.F90 :
344
345	IF( kt == nit000) THEN
346	!--- computation of rdtbs2 ===> !gm is it really usefull ????
347	IF( nacc == 1 ) THEN
348	rdtbs2 = nfbulk * rdtmin * 0.5
349	ELSE
350	rdtbs2 = nfbulk * rdt * 0.5
351	ENDIF
352	IF( .NOT.ln_rstart ) THEN
353	zcof = REAL( nfbulk - 1, wp )
354	gsst(:,:) = zcof * ( tn(:,:,1) + rt0 )
355	gsss(:,:) = zcof * tn_ice(:,:)
356	gu (:,:) = zcof * un(:,:,1)
357	gv (:,:) = zcof * vn(:,:,1)
358	ENDIF
359	ENDIF
360
361	! ----------------------------------------------------------------------------- !
362	! 0 Mean SST, ice temperature and ocean currents over nfbulk pdt !
363	! ----------------------------------------------------------------------------- !
364	!
365	gsst(:,:) = gsst(:,:) + (tn(:,:,1) + rt0 )
366	gsss(:,:) = gsss(:,:) + tn_ice(:,:)
367	gu (:,:) = gu (:,:) + un (:,:,1)
368	gv (:,:) = gv (:,:) + vn (:,:,1)
369
370	IF( MOD( kt - 1 , nfbulk ) == 0 ) THEN
371	!
372	zsst(:,:) = gsst(:,:) / REAL( nfbulk, wp ) * tmask(:,:,1) ! mean sst in K
373	zsss(:,:) = gsss(:,:) / REAL( nfbulk ) * tmask(:,:,1) ! mean tn_ice in K
374
375	where ( zsst(:,:) .eq. 0.e0 ) zsst(:,:) = rt0 !lb vilain !!???
376	where ( zsss(:,:) .eq. 0.e0 ) zsss(:,:) = rt0 !lb //
377
378	!!gm better coded:
379	! Caution set to rt0 over land, not 0 Kelvin (otherwise floating point exception in bulk computation)
380	! zcof = 1. / REAL( nfbulk, wp )
381	! zsst(:,:) = gsst(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1) ! mean sst in K
382	! zsss(:,:) = gsss(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1) ! mean tn_ice in K
383	!!gm
384
385	zut(:,:) = 0.e0 !!gm not necessary but at least first and last column
386	zvt(:,:) = 0.e0
387	! lb
388	! Interpolation of surface current at T-point, zut and zvt :
389	DO ji=2, jpi
390	zut(ji,:) = 0.5*(gu(ji-1,:) + gu(ji,:)) / REAL( nfbulk )
391	END DO
392	!
393	DO jj=2, jpj
394	zvt(:,jj) = 0.5*(gv(:,jj-1) + gv(:,jj)) / REAL( nfbulk )
395	END DO
396	!!gm better coded
397	! zcof = 0.5 / REAL( nfbulk, wp )
398	! DO jj = 2, jpjm1
399	! DO ji = fs_2, fs_jpim1 ! vector opt.
400	! zut(ji,jj) = ( gu(ji-1,jj ) + gu(ji,jj) ) * zcof
401	! zvt(ji,jj) = ( gv(ji ,jj-1) + gv(ji,jj) ) * zcof
402	! END DO
403	! END DO
404	!!gm
405
406	CALL lbc_lnk( zut(:,:), 'T', -1. )
407	CALL lbc_lnk( zvt(:,:), 'T', -1. )
408
409	! ----------------------------------------------------------------------------- !
410	! I Radiative FLUXES !
411	! ----------------------------------------------------------------------------- !
412
413	!--- Ocean & Ice Albedo
414	alb_oce_os(:,:) = 0. ; alb_oce_cs(:,:) = 0. !gm is it necessary ???
415	alb_ice_os(:,:) = 0. ; alb_ice_cs(:,:) = 0.
416	CALL flx_blk_albedo( alb_ice_os, alb_oce_os, alb_ice_cs, alb_oce_cs )
417
418	!--- Radiative fluxes over ocean
419	qsr_oce(:,:) = ( 1. - 0.066 ) * flxnow(:,:,5) ! Solar Radiation
420	qlw_oce(:,:) = flxnow(:,:,6) - Stefzsst(:,:)zsst(:,:)zsst(:,:)zsst(:,:) ! Long Waves (Infra-red)
421
422	!--- Radiative fluxes over ice
423	qsr_ice(:,:) = ( 1. - alb_ice_cs(:,:) ) * flxnow(:,:,5) ! Solar Radiation
424	qlw_ice(:,:) = 0.95 * ( flxnow(:,:,6) - Stefzsss(:,:)zsss(:,:)zsss(:,:)zsss(:,:) ) ! Long Waves (Infra-red)
425	!-------------------------------------------------------------------------------!
426
427
428	! ----------------------------------------------------------------------------- !
429	! II Turbulent FLUXES !
430	! ----------------------------------------------------------------------------- !
431	!
432	! scalar wind ( = \| U10m - SSU \| )
433	! It is important to take into account the sea surface courant
434	! lb
435	! Now, wind components are provided on T-points within the netcdf input file.
436	! Thus, the wind module is computded at T-points taking into account the sea
437
438	!--- Module of relative wind
439	dUnormt(:,:) = sqrt( (flxnow(:,:,2) - zut(:,:))*(flxnow(:,:,2) - zut(:,:)) &
440	& + (flxnow(:,:,3) - zvt(:,:))*(flxnow(:,:,3) - zvt(:,:)) )
441	!!gm more efficient coding:
442	! DO jj = 1, jpi
443	! DO ji = 1, jpi
444	! zzu = flxnow(ji,jj,2) - zut(ji,jj)
445	! zzv = flxnow(ji,jj,3) - zvt(ji,jj)
446	! dUnormt(ji,jj) = SQRT( zzuzzu + zzvzzv )
447	! END DO
448	! END DO
449	!!gm
450	! lb
451	!
452	!--- specific humidity at temp SST
453	qsatw(:,:) = 0.98640380exp(-5107.4/zsst(:,:))/rhoa
454	!
455	!--- specific humidity at temp tn_ice
456	qsat(:,:) = 11637800*exp(-5897.8/zsss(:,:))/rhoa
457	!
458
459	! CORE iterartive algo for computation of Cd, Ch, Ce at T-point :
460	! ===============================================================
461	IF( ln_2m ) THEN
462	IF( kt == nit000 .AND. lwp ) THEN
463	WRITE(numout,*)
464	WRITE(numout,*)' flx_core : global CORE fields in NetCDF format'
465	WRITE(numout,*)' ~~~~~~~~~ '
466	WRITE(numout,*) ' Calling TURB_CORE_2Z for bulk transfert coefficients'
467	WRITE(numout,*)
468	ENDIF
469	!! If air temp. and spec. hum. are given at different height (2m) than wind (10m) :
470	CALL TURB_CORE_2Z(2.,10.,zsst(:,:), flxnow(:,:,7), qsatw(:,:), flxnow(:,:,4), &
471	& dUnormt(:,:), Cd(:,:), Ch(:,:), Ce(:,:), zt_zu(:,:), zq_zu(:,:))
472	ELSE
473	IF( kt == nit000 .AND. lwp ) THEN
474	WRITE(numout,*)
475	WRITE(numout,*)' flx_core : global CORE fields in NetCDF format'
476	WRITE(numout,*)' ~~~~~~~~~~ '
477	WRITE(numout,*) ' Calling TURB_CORE_1Z for bulk transfert coefficients'
478	WRITE(numout,*)
479	ENDIF
480	!! If air temp. and spec. hum. are given at same height than wind (10m) :
481	CALL TURB_CORE_1Z(10., zsst(:,:), flxnow(:,:,7), qsatw(:,:), flxnow(:,:,4), &
482	& dUnormt(:,:), Cd(:,:), Ch(:,:), Ce(:,:) )
483	ENDIF
484
485	! II.1) Momentum over ocean and ice
486	! ---------------------------------
487	! Tau_x at T-point
488	tauxt(:,:) = rhoadUnormt(:,:)( (1. - freeze(:,:))Cd(:,:)(flxnow(:,:,2) - zut(:,:)) &
489	& + freeze(:,:)Ciceflxnow(:,:,2) ) !lb correct pour glace
490	! Tau_y at T-point
491	tauyt(:,:) = rhoadUnormt(:,:)( (1. - freeze(:,:))Cd(:,:)(flxnow(:,:,3) - zvt(:,:)) &
492	& + freeze(:,:)Ciceflxnow(:,:,3) )
493	!
494	CALL lbc_lnk( tauxt(:,:), 'T', -1. ) !!gm seems not decessary at this point....
495	CALL lbc_lnk( tauyt(:,:), 'T', -1. )
496	!
497
498	! Katabatic winds parametrization
499	! -------------------------------
500	IF( ln_kata ) THEN
501	!
502	IF( kt == nit000 ) THEN
503	IF (lwp ) WRITE(numout,*) ' Katabatic winds parametrization is active'
504	clnamekata = 'katamask.nc'
505	clvarkatax = 'katamaskx'
506	clvarkatay = 'katamasky'
507
508	#if defined key_agrif
509	IF ( .NOT. Agrif_Root() ) THEN
510	clnamekata = TRIM(Agrif_CFixed())//'_'//TRIM(clnamekata)
511	ENDIF
512	#endif
513	CALL iom_open( TRIM(clnamekata), inum )
514	CALL iom_get ( inum, jpdom_data, TRIM(clvarkatax), rmskkatax )
515	CALL iom_get ( inum, jpdom_data, TRIM(clvarkatay), rmskkatay )
516	CALL iom_close ( inum )
517
518	WHERE( (rmskkatax < 0.000001) .AND. (rmskkatax > -0.000001) )
519	rmskkatax=1
520	rmskkatay=1
521	END WHERE
522
523	CALL lbc_lnk( rmskkatax(:,:), 'T', 1. ) ; CALL lbc_lnk( rmskkatay(:,:), 'T', 1. )
524
525	IF(MAXVAL(rmskkatax) > 6.00001) THEN
526	WRITE(ctmp1,*) 'Problem in the data mask : maxval = ',MAXVAL(rmskkatax),' ( it is GT 6)'
527	CALL ctl_stop( ctmp1 )
528	ENDIF
529	ENDIF
530
531	! apply katabatic wind enhancement on grid T (before projection)
532	tauxt(:,:) = rmskkatax(:,:) * tauxt(:,:)
533	tauyt(:,:) = rmskkatay(:,:) * tauyt(:,:)
534	!
535	ENDIF
536	!
537	!Tau_x at U-point
538	DO ji = 1, jpim1
539	taux(ji,:) = 0.5*(tauxt(ji,:) + tauxt(ji+1,:))
540	END DO
541	!
542	! Tau_y at V-point
543	DO jj = 1, jpjm1
544	tauy(:,jj) = 0.5*(tauyt(:,jj) + tauyt(:,jj+1))
545	END DO
546
547	!!gm ==> at this point, a lbc is required, no?
548	!
549	! lb : should we do this here? !!gm yes should we remove that???
550	tauxg(:,:) = taux(:,:) ! Save components in
551	tauyg(:,:) = tauy(:,:) ! geographical ref on U grid
552	!
553	!
554	! II.2) Turbulent fluxes over ocean
555	! ---------------------------------
556	!
557	IF ( ln_2m ) THEN
558	!!
559	!! Values of temp. and hum. adjusted to 10m must be used instead of 2m values
560	!! Sensible Heat : ! right sign for ocean
561	qsb_oce(:,:) = rhoacpCh(:,:)(zt_zu(:,:) - zsst(:,:))dUnormt(:,:)
562	!!
563	!! Latent Heat : ! wrong sign for ocean
564	evap(:,:) = rhoaCe(:,:)(qsatw(:,:) - zq_zu(:,:))*dUnormt(:,:)
565	!!
566	ELSE
567	!!
568	!! Sensible Heat : ! right sign for ocean
569	qsb_oce(:,:) = rhoacpCh(:,:)(flxnow(:,:,7) - zsst(:,:))dUnormt(:,:)
570	!!
571	!! Latent Heat : ! wrong sign for ocean
572	evap(:,:) = rhoaCe(:,:)(qsatw(:,:) - flxnow(:,:,4))*dUnormt(:,:)
573	!!
574	END IF
575	!!
576	qla_oce(:,:) = Lv * evap(:,:) ! right sign for ocean
577
578	!
579	! II.3) Turbulent fluxes over ice
580	! -------------------------------
581	!
582	! Sensible Heat :
583	qsb_ice(:,:) = rhoacpCice(zsss(:,:) - flxnow(:,:,7))dUnormt(:,:) !lb use
584	! !lb dUnormt???
585	! Latent Heat : !lb or rather Unormt?
586	qla_ice(:,:) = LsrhoaCice(qsat(:,:) - flxnow(:,:,4))dUnormt(:,:)
587	!
588	!-------------------------------------------------------------------------------------
589
590
591
592	! ----------------------------------------------------------------------------- !
593	! III Total FLUXES !
594	! ----------------------------------------------------------------------------- !
595	!
596	! III.1) Downward Non Solar flux over ocean
597	! -----------------------------------------
598	qnsr_oce(:,:) = qlw_oce(:,:) - qsb_oce(:,:) - qla_oce(:,:)
599	!
600	! III.1) Downward Non Solar flux over ice
601	! ---------------------------------------
602	qnsr_ice(:,:) = qlw_ice(:,:) - qsb_ice(:,:) - qla_ice(:,:)
603	!
604	!-------------------------------------------------------------------------------!
605
606
607
608	! 6. TOTAL NON SOLAR SENSITIVITY
609
610	dqlw_ice(:,:)= 4.00.95Stefzsss(:,:)zsss(:,:)*zsss(:,:)
611
612	! d qla_ice/ d zsss
613	dqla_ice(:,:) = -LsCice0.9811637800/(rhoarhoa) &
614	* (-5897.8)/(zsss(:,:)*zsss(:,:)) &
615	* exp(-5897.8/zsss(:,:)) * dUnormt(:,:)
616
617	! d qsb_ice/ d zsss
618	dqsb_ice(:,:) = rhoa * cp * Cice * dUnormt(:,:)
619
620	dqns_ice(:,:) = - ( dqlw_ice(:,:) + dqsb_ice(:,:) + dqla_ice(:,:) )
621
622
623	!--------------------------------------------------------------------
624	! FRACTION of net shortwave radiation which is not absorbed in the
625	! thin surface layer and penetrates inside the ice cover
626	! ( Maykut and Untersteiner, 1971 ; Elbert and Curry, 1993 )
627
628	catm1(:,:) = 1.0 - 0.3 ! flxnow(:,:,8)
629	fr1_i0(:,:) = &
630	(0.18 * catm1(:,:) + 0.35 * flxnow(:,:,8))
631	fr2_i0(:,:) = &
632	(0.82 * catm1(:,:) + 0.65 * flxnow(:,:,8))
633
634
635	! Total PRECIPITATION (kg/m2/s)
636	tprecip(:,:) = flxnow(:,:,1)
637	! rename precipitation for freshwater budget calculations:
638	watm(:,:) = flxnow(:,:,1)*1000
639	!
640
641	! SNOW PRECIPITATION (kg/m2/s)
642	sprecip(:,:) = flxnow(:,:,8)
643
644	!---------------------------------------------------------------------
645
646
647
648	CALL lbc_lnk( taux (:,:) , 'U', -1. )
649	CALL lbc_lnk( tauy (:,:) , 'V', -1. )
650	CALL lbc_lnk( qsr_oce (:,:) , 'T', 1. )
651	CALL lbc_lnk( qnsr_oce(:,:) , 'T', 1. )
652	CALL lbc_lnk( qsr_ice (:,:) , 'T', 1. )
653	CALL lbc_lnk( qnsr_ice(:,:) , 'T', 1. )
654	CALL lbc_lnk( qla_ice (:,:) , 'T', 1. )
655	CALL lbc_lnk( dqns_ice(:,:) , 'T', 1. )
656	CALL lbc_lnk( dqla_ice(:,:) , 'T', 1. )
657	CALL lbc_lnk( fr1_i0 (:,:) , 'T', 1. )
658	CALL lbc_lnk( fr2_i0 (:,:) , 'T', 1. )
659	CALL lbc_lnk( tprecip (:,:) , 'T', 1. )
660	CALL lbc_lnk( sprecip (:,:) , 'T', 1. )
661	CALL lbc_lnk( evap (:,:) , 'T', 1. )
662	!!
663	!!
664	!! NEVER mask the windstress!!
665	qsr_oce (:,:) = qsr_oce (:,:)*tmask(:,:,1)
666	qnsr_oce(:,:) = qnsr_oce(:,:)*tmask(:,:,1)
667	qsr_ice (:,:) = qsr_ice (:,:)*tmask(:,:,1)
668	qnsr_ice(:,:) = qnsr_ice(:,:)*tmask(:,:,1)
669	qla_ice (:,:) = qla_ice (:,:)*tmask(:,:,1)
670	dqns_ice(:,:) = dqns_ice(:,:)*tmask(:,:,1)
671	dqla_ice(:,:) = dqla_ice(:,:)*tmask(:,:,1)
672	fr1_i0 (:,:) = fr1_i0 (:,:)*tmask(:,:,1)
673	fr2_i0 (:,:) = fr2_i0 (:,:)*tmask(:,:,1)
674	tprecip (:,:) = tprecip (:,:)*tmask(:,:,1)
675	sprecip (:,:) = sprecip (:,:)*tmask(:,:,1)
676	evap (:,:) = evap (:,:)*tmask(:,:,1)
677	gsst(:,:) = 0.e0
678	gsss(:,:) = 0.e0
679	gu (:,:) = 0.e0
680	gv (:,:) = 0.e0
681
682	! Degugging print (A.M. Treguier)
683	! write (55) taux, tauy, qnsr_oce, qsb_ice, qnsr_ice, qla_ice, dqns_ice &
684	! & , dqla_ice, fr1_i0, fr2_i0, qlw_oce, qla_oce, qsb_oce, evap
685	! write(numout,*) 'written 14 arrays on unit fort.55'
686
687
688	ENDIF
689
690	! ------------------- !
691	! Last call kt=nitend !
692	! ------------------- !
693
694	! Closing of the numflx file (required in mpp)
695
696	IF( kt == nitend ) THEN
697	DO ji=1, jpfile
698	CALL iom_close(numflxall(ji))
699	END DO
700	ENDIF
701
702	END SUBROUTINE flx
703
704
705	SUBROUTINE TURB_CORE_1Z( zzu, T_0, T_a, q_sat, q_a, &
706	& dU , C_d, C_h , C_e )
707	!!----------------------------------------------------------------------
708	!! * ROUTINE turb_core *
709	!!
710	!! ** Purpose : Computes turbulent transfert coefficients of surface
711	!! fluxes according to Large & Yeager (2004).
712	!!
713	!! ** 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
714	!! Momentum, Latent and sensible heat exchange coefficients
715	!! Caution: this procedure should only be used in cases when air
716	!! temperature (T_air), air specific humidity (q_air) and wind (dU)
717	!! are provided at the same height 'zzu'!
718	!!
719	!! References :
720	!! Large & Yeager, 2004 : ???
721	!! History :
722	!! ! XX-XX (??? ) Original code
723	!! 9.0 ! 05-08 (L. Brodeau) Optimisation
724	!!----------------------------------------------------------------------
725	!! * Arguments
726	REAL(wp), INTENT( in ) :: &
727	zzu ! height for all given atmospheric variables [m]
728	REAL(wp), INTENT( in ), DIMENSION(jpi,jpj) :: &
729	T_0, & ! sea surface temperature [Kelvin]
730	T_a, & ! potential air temperature at zzu [Kelvin]
731	q_sat, & ! sea surface specific humidity at zzu [kg/kg]
732	q_a, & ! specific air humidity at zzu [kg/kg]
733	dU ! relative wind module \|U(zzu)-U(0)\| at zzu [m/s]
734	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: &
735	C_d, & ! transfer coefficient for momentum (tau)
736	C_h, & ! transfer coefficient for sensible heat (Q_sens)
737	C_e ! tansfert coefficient for evaporation (Q_lat)
738
739	!! * Local declarations
740	REAL(wp), DIMENSION(jpi,jpj) :: &
741	dU10, & ! dU [m/s]
742	dT, & ! air/sea temperature differeence [K]
743	dq, & ! air/sea humidity difference [K]
744	Cd_n10, & ! 10m neutral drag coefficient
745	Ce_n10, & ! 10m neutral latent coefficient
746	Ch_n10, & ! 10m neutral sensible coefficient
747	Cd, & ! drag coefficient
748	Ce, & ! latent coefficient
749	Ch, & ! sensible coefficient
750	sqrt_Cd_n10, & ! root square of Cd_n10
751	sqrt_Cd, & ! root square of Cd
752	T_vpot, & ! virtual potential temperature [K]
753	T_star, & ! turbulent scale of tem. fluct.
754	q_star, & ! turbulent humidity of temp. fluct.
755	U_star, & ! turb. scale of velocity fluct.
756	L, & ! Monin-Obukov length [m]
757	zeta, & ! stability parameter at height zzu
758	X2, X, &
759	psi_m, &
760	psi_h, &
761	U_n10, & ! neutral wind velocity at 10m [m]
762	xlogt
763
764	INTEGER :: jk ! doomy loop for iterations
765	INTEGER, PARAMETER :: nit = 3 ! number of iterations
766
767	INTEGER, DIMENSION(jpi,jpj) :: &
768	stab, & ! stability test integer
769	stabit ! stability within iterative loop
770
771	REAL(wp), PARAMETER :: &
772	pi = 3.14159, & ! Pi
773	grav = 9.8, & ! gravity
774	kappa = 0.4 ! von Karman's constant
775	!!----------------------------------------------------------------------
776	!!
777	!! -------------
778	!! S T A R T
779	!! -------------
780	!!
781	!! I.1 ) Preliminary stuffs
782	!! ------------------------
783	!!
784	!! Air/sea differences
785	!! -------------------
786	dU10 = max(0.5, dU) ! we don't want to fall under 0.5 m/s
787	dT = T_a - T_0 ! assuming that T_a is allready the potential temp. at zzu
788	dq = q_a - q_sat
789	!!
790	!! Virtual potential temperature
791	!! -----------------------------
792	T_vpot = T_a(1. + 0.608q_a)
793	!!
794	!!
795	!! I.2 ) Computing Neutral Drag Coefficient
796	!! ----------------------------------------
797	Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! \\ L & Y eq. (6a)
798	sqrt_Cd_n10 = sqrt(Cd_n10)
799	!!
800	Ce_n10 = 1E-3 * ( 34.6 * sqrt_Cd_n10 ) ! \\ L & Y eq. (6b)
801	!!
802	!! First guess of stabilitty :
803	stab = 0.5 + sign(0.5,dT) ! stable : stab = 1 ; unstable : stab = 0
804	Ch_n10 = 1E-3sqrt_Cd_n10(18stab + 32.7(1-stab)) ! \\ L & Y eq. (6c), (6d)
805	!!
806	!! Initializing transfert coefficients with their first guess neutral equivalents :
807	Cd = Cd_n10 ; Ce = Ce_n10 ; Ch = Ch_n10
808	!!
809	!!
810	!!
811	!! II ) Now starting iteration loop (IDM)
812	!! ----------------------------------------
813	!!
814	DO jk=1, nit
815	!!
816	sqrt_Cd = sqrt(Cd)
817	!!
818	!! Turbulent scales :
819	!! ------------------
820	U_star = sqrt_Cd*dU10 ! \\ L & Y eq. (7a)
821	T_star = Ch/sqrt_Cd*dT ! \\ L & Y eq. (7b)
822	q_star = Ce/sqrt_Cd*dq ! \\ L & Y eq. (7c)
823	!!
824	!! Estimate the Monin-Obukov length :
825	!! ----------------------------------
826	!dbg1 = T_star/T_vpot
827	!dbg2 = q_star/(q_a + 1./0.608)
828	L = (U_star*2)/( kappagrav*(T_star/T_vpot + q_star/(q_a + 1./0.608)) )
829	!!
830	!! Stability parameters :
831	!! ----------------------
832	zeta = zzu/L
833	zeta = sign( min(abs(zeta),10.0), zeta )
834	!!
835	!! Psis, L & Y eq. (8c), (8d), (8e) :
836	!! ----------------------------------
837	X2 = sqrt(abs(1. - 16.*zeta)) ; X2 = max(X2 , 1.0) ; X = sqrt(X2)
838	!!
839	stabit = 0.5 + sign(0.5,zeta)
840	!!
841	! dbg1 = -5zetastabit
842	! dbg2 = 2*log((1. + X)/2)
843	! dbg3 = log((1. + X2)/2)
844	! dbg4 = 2*atan(X)
845
846	psi_m = -5zetastabit & ! Stable
847	& + (1 - stabit)(2log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable
848	!!
849	psi_h = -5zetastabit & ! Stable
850	& + (1 - stabit)(2log( (1. + X2)/2 )) ! Unstable
851	!!
852	!!
853	!! Shifting the wind speed to 10m and neutral stability :
854	!! ------------------------------------------------------
855	U_n10 = dU101./(1. + sqrt(Cd_n10)/kappa(log(zzu/10.) - psi_m)) ! \\ L & Y eq. (9a)
856	!!
857	!!
858	!! Updating the neutral 10m transfer coefficients :
859	!! ------------------------------------------------
860	Cd_n10 = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! \\ L & Y eq. (6a)
861	sqrt_Cd_n10 = sqrt(Cd_n10)
862	!!
863	Ce_n10 = 1E-3 * (34.6 * sqrt_Cd_n10) ! \\ L & Y eq. (6b)
864	!!
865	stab = 0.5 + sign(0.5,zeta)
866	Ch_n10 = 1E-3sqrt_Cd_n10(18.stab + 32.7(1-stab)) ! \\ L & Y eq. (6c), (6d)
867	!!
868	!!
869	!! Shifting the neutral 10m transfer coefficients to ( zzu , zeta ) :
870	!! --------------------------------------------------------------------
871	!! Problem here, formulation used within L & Y differs from the one provided
872	!! in their fortran code (only for Ce and Ch)
873	!!
874	Cd = Cd_n10/(1. + sqrt_Cd_n10/kappa(log(zzu/10) - psi_m))*2 ! \\ L & Y eq. (10a)
875	!!
876	xlogt = log(zzu/10) - psi_h
877	!? Ch = Ch_n10sqrt(Cd/Cd_n10)/(1. + Ch_n10/(kappasqrt_Cd_n10)*xlogt)
878	Ch = Ch_n10/( 1. + Ch_n10xlogt/kappa/sqrt_Cd_n10 )*2 ! \\ L & Y eq. (10b)
879	!!
880	!? Ce = Ce_n10sqrt(Cd/Cd_n10)/(1. + Ce_n10/(kappasqrt_Cd_n10)*xlogt)
881	Ce = Ce_n10/( 1. + Ce_n10xlogt/kappa/sqrt_Cd_n10 )*2 ! \\ L & Y eq. (10c)
882	!!
883	!!
884	END DO
885	!!
886	C_d = Cd
887	C_h = Ch
888	C_e = Ce
889	!!
890	END SUBROUTINE TURB_CORE_1Z
891
892
893	SUBROUTINE TURB_CORE_2Z(zt, zu, sst, T_zt, q_sat, q_zt, dU, Cd, Ch, Ce, T_zu, q_zu)
894	!!----------------------------------------------------------------------
895	!! * ROUTINE turb_core *
896	!!
897	!! ** Purpose : Computes turbulent transfert coefficients of surface
898	!! fluxes according to Large & Yeager (2004).
899	!!
900	!! ** 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
901	!! Momentum, Latent and sensible heat exchange coefficients
902	!! Caution: this procedure should only be used in cases when air
903	!! temperature (T_air) and air specific humidity (q_air) are at 2m
904	!! whereas wind (dU) is at 10m.
905	!!
906	!! References :
907	!! Large & Yeager, 2004 : ???
908	!! History :
909	!! ! XX-XX (??? ) Original code
910	!! 9.0 ! 05-08 (L. Brodeau) Optimisation
911	!!----------------------------------------------------------------------
912	!! * Arguments
913	REAL(wp), INTENT( in ) :: &
914	zt, & ! height for T_zt and q_zt [m]
915	zu ! height for dU [m]
916	!!
917	REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: &
918	sst, & ! sea surface temperature [Kelvin]
919	T_zt, & ! potential air temperature [Kelvin]
920	q_sat, & ! sea surface specific humidity [kg/kg]
921	q_zt, & ! specific air humidity [kg/kg]
922	dU ! wind module \|U(zu)-U(0)\| [m/s]
923	!!
924	REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: &
925	Cd, Ch, Ce, &
926	T_zu, & ! air temp. shifted at zu [K]
927	q_zu ! spec. hum. shifted at zu [kg/kg]
928
929	!! * Local declarations
930	REAL(wp), DIMENSION(jpi,jpj) :: &
931	dU10, & ! dU [m/s]
932	dT, & ! air/sea temperature differeence [K]
933	dq, & ! air/sea humidity difference [K]
934	Cd_n10, & ! 10m neutral drag coefficient
935	Ce_n10, & ! 10m neutral latent coefficient
936	Ch_n10, & ! 10m neutral sensible coefficient
937	sqrt_Cd_n10, & ! root square of Cd_n10
938	sqrt_Cd, & ! root square of Cd
939	T_vpot_u, & ! virtual potential temperature [K]
940	T_star, & ! turbulent scale of tem. fluct.
941	q_star, & ! turbulent humidity of temp. fluct.
942	U_star, & ! turb. scale of velocity fluct.
943	L, & ! Monin-Obukov length [m]
944	zeta_u, & ! stability parameter at height zu
945	zeta_t, & ! stability parameter at height zt
946	U_n10, & ! neutral wind velocity at 10m [m]
947	xlogt, xct
948	!!
949	INTEGER, PARAMETER :: nb_itt = 3 ! number of itterations
950	!!
951	!! Some physical constants :
952	REAL(wp), PARAMETER :: &
953	grav = 9.8, & ! gravity
954	kappa = 0.4 ! von Karman's constant
955	!!
956	INTEGER :: j_itt
957	INTEGER, DIMENSION(jpi,jpj) :: &
958	stab ! 1st stability test integer
959	!!----------------------------------------------------------------------
960	!!
961	!! -------------
962	!! S T A R T
963	!! -------------
964	!!
965	!! I.1 ) Preliminary stuffs
966	!! ------------------------
967	!!
968	!! Initial air/sea differences
969	dU10 = max(0.5, dU) ; dT = T_zt - sst ; dq = q_zt - q_sat
970	!!
971	!! Neutral Drag Coefficient :
972	stab = 0.5 + sign(0.5,dT) ! stab = 1 if dT > 0 -> STABLE
973	Cd_n10 = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 )
974	sqrt_Cd_n10 = sqrt(Cd_n10)
975	Ce_n10 = 1E-3( 34.6 sqrt_Cd_n10 )
976	Ch_n10 = 1E-3sqrt_Cd_n10(18stab + 32.7(1 - stab))
977	!!
978	!! I.2 ) Computing Neutral Drag Coefficient
979	!! ----------------------------------------
980	!! Initializing transf. coeff. with their first guess neutral equivalents :
981	Cd = Cd_n10 ; Ce = Ce_n10 ; Ch = Ch_n10 ; sqrt_Cd = sqrt(Cd)
982	!!
983	!! Initializing z_u values with z_t values :
984	T_zu = T_zt ; q_zu = q_zt
985	!!
986	!! II ) Now starting iteration loop (IDM)
987	!! ----------------------------------------
988	!!
989	DO j_itt=1, nb_itt
990	!!
991	!! Updating air/sea differences :
992	dT = T_zu - sst ; dq = q_zu - q_sat
993	!!
994	!! Updating virtual potential temperature at zu :
995	T_vpot_u = T_zu(1. + 0.608q_zu)
996	!!
997	!! Updating turbulent scales : (L & Y eq. (7))
998	U_star = sqrt_CddU10 ; T_star = Ch/sqrt_CddT ; q_star = Ce/sqrt_Cd*dq
999	!!
1000	!! Estimate the Monin-Obukov length at height zu :
1001	L = (U_star*U_star) &
1002	& / (kappagrav/T_vpot_u(T_star(1.+0.608q_zu) + 0.608T_zuq_star))
1003	!!
1004	!! Stability parameters :
1005	zeta_u = zu/L ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u )
1006	zeta_t = zt/L ; zeta_t = sign( min(abs(zeta_t),10.0), zeta_t )
1007	!!
1008	!! Shifting the wind speed to 10m and neutral stability : (L & Y eq.(9a))
1009	U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u)))
1010	!!
1011	!! Shifting temperature and humidity at zu : (L & Y eq. (9b-9c))
1012	T_zu = T_zt - T_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t))
1013	q_zu = q_zt - q_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t))
1014	!!
1015	!! q_zu cannot have a negative value : forcing 0
1016	stab = 0.5 + sign(0.5,q_zu) ; q_zu = stab*q_zu
1017	!!
1018	!! Updating the neutral 10m transfer coefficients :
1019	Cd_n10 = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a)
1020	sqrt_Cd_n10 = sqrt(Cd_n10)
1021	Ce_n10 = 1E-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b)
1022	stab = 0.5 + sign(0.5,zeta_u)
1023	Ch_n10 = 1E-3sqrt_Cd_n10(18.stab + 32.7(1-stab)) ! L & Y eq. (6c-6d)
1024	!!
1025	!!
1026	!! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) :
1027	xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u))
1028	Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd)
1029	!!
1030	xlogt = log(zu/10.) - psi_h(zeta_u)
1031	!!
1032	xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10
1033	Ch = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct
1034	!!
1035	xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10
1036	Ce = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct
1037	!!
1038	!!
1039	END DO
1040	!!
1041	END SUBROUTINE TURB_CORE_2Z
1042
1043	FUNCTION psi_m(zta) !! Psis, L & Y eq. (8c), (8d), (8e)
1044	REAL(wp), PARAMETER :: pi = 3.14159
1045	REAL(wp), DIMENSION(jpi,jpj) :: psi_m
1046	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta
1047	REAL(wp), DIMENSION(jpi,jpj) :: X2, X, stabit
1048	X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.0) ; X = sqrt(X2)
1049	stabit = 0.5 + sign(0.5,zta)
1050	psi_m = -5.ztastabit & ! Stable
1051	& + (1. - stabit)(2log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable
1052	END FUNCTION psi_m
1053
1054	FUNCTION psi_h(zta) !! Psis, L & Y eq. (8c), (8d), (8e)
1055	REAL(wp), DIMENSION(jpi,jpj) :: psi_h
1056	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta
1057	REAL(wp), DIMENSION(jpi,jpj) :: X2, X, stabit
1058	X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.) ; X = sqrt(X2)
1059	stabit = 0.5 + sign(0.5,zta)
1060	psi_h = -5.ztastabit & ! Stable
1061	& + (1. - stabit)(2.log( (1. + X2)/2. )) ! Unstable
1062	END FUNCTION psi_h
1063
1064
1065	SUBROUTINE flx_blk_albedo( palb , palcn , palbp , palcnp )
1066	!!----------------------------------------------------------------------
1067	!! * ROUTINE flx_blk_albedo *
1068	!!
1069	!! ** Purpose : Computation of the albedo of the snow/ice system
1070	!! as well as the ocean one
1071	!!
1072	!! ** Method : - Computation of the albedo of snow or ice (choose the
1073	!! right one by a large number of tests
1074	!! - Computation of the albedo of the ocean
1075	!!
1076	!! References :
1077	!! Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250.
1078	!!
1079	!! History :
1080	!! 8.0 ! 01-04 (LIM 1.0)
1081	!! 8.5 ! 03-07 (C. Ethe, G. Madec) Optimization (old name:shine)
1082	!!----------------------------------------------------------------------
1083	!! * Modules used
1084	USE ice ! ???
1085
1086	!! * Arguments
1087	REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: &
1088	palb , & ! albedo of ice under overcast sky
1089	palcn , & ! albedo of ocean under overcast sky
1090	palbp , & ! albedo of ice under clear sky
1091	palcnp ! albedo of ocean under clear sky
1092
1093	!! * Local variables
1094	INTEGER :: &
1095	ji, jj ! dummy loop indices
1096	REAL(wp) :: &
1097	c1 = 0.05 , & ! constants values
1098	c2 = 0.1 , &
1099	albice = 0.50 , & ! albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
1100	cgren = 0.06 , & ! correction of the snow or ice albedo to take into account
1101	! effects of cloudiness (Grenfell & Perovich, 1984)
1102	alphd = 0.80 , & ! coefficients for linear interpolation used to compute
1103	alphdi = 0.72 , & ! albedo between two extremes values (Pyane, 1972)
1104	alphc = 0.65 , &
1105	zmue = 0.4 , & ! cosine of local solar altitude
1106	zzero = 0.0 , &
1107	zone = 1.0
1108
1109	REAL(wp) :: &
1110	zmue14 , & ! zmue**1.4
1111	zalbpsnm , & ! albedo of ice under clear sky when snow is melting
1112	zalbpsnf , & ! albedo of ice under clear sky when snow is freezing
1113	zalbpsn , & ! albedo of snow/ice system when ice is coverd by snow
1114	zalbpic , & ! albedo of snow/ice system when ice is free of snow
1115	zithsn , & ! = 1 for hsn >= 0 ( ice is cov. by snow ) ; = 0 otherwise (ice is free of snow)
1116	zitmlsn , & ! = 1 freezinz snow (sist >=rt0_snow) ; = 0 melting snow (sist<rt0_snow)
1117	zihsc1 , & ! = 1 hsn <= c1 ; = 0 hsn > c1
1118	zihsc2 ! = 1 hsn >= c2 ; = 0 hsn < c2
1119	REAL(wp), DIMENSION(jpi,jpj) :: &
1120	zalbfz , & ! ( = alphdi for freezing ice ; = albice for melting ice )
1121	zficeth ! function of ice thickness
1122	LOGICAL , DIMENSION(jpi,jpj) :: &
1123	llmask
1124	!! to be included for without seaice
1125	!! REAL(wp), DIMENSION(jpi,jpj) :: & !:
1126	!! sist , & !: Sea-Ice Surface Temperature (Kelvin )
1127	!! hsnif , & !: Snow thickness
1128	!! hicif !: Ice thickness
1129
1130	!!---------------------------------------------------------------------
1131
1132	!-------------------------
1133	! Computation of zficeth
1134	!--------------------------
1135
1136	llmask = (hsnif == 0.0) .AND. ( sist >= rt0_ice )
1137	WHERE ( llmask ) ! ice free of snow and melts
1138	zalbfz = albice
1139	ELSEWHERE
1140	zalbfz = alphdi
1141	END WHERE
1142
1143	DO jj = 1, jpj
1144	DO ji = 1, jpi
1145	IF( hicif(ji,jj) > 1.5 ) THEN
1146	zficeth(ji,jj) = zalbfz(ji,jj)
1147	ELSEIF( hicif(ji,jj) > 1.0 .AND. hicif(ji,jj) <= 1.5 ) THEN
1148	zficeth(ji,jj) = 0.472 + 2.0 * ( zalbfz(ji,jj) - 0.472 ) * ( hicif(ji,jj) - 1.0 )
1149	ELSEIF( hicif(ji,jj) > 0.05 .AND. hicif(ji,jj) <= 1.0 ) THEN
1150	zficeth(ji,jj) = 0.2467 + 0.7049 * hicif(ji,jj) &
1151	& - 0.8608 * hicif(ji,jj) * hicif(ji,jj) &
1152	& + 0.3812 * hicif(ji,jj) * hicif(ji,jj) * hicif (ji,jj)
1153	ELSE
1154	zficeth(ji,jj) = 0.1 + 3.6 * hicif(ji,jj)
1155	ENDIF
1156	END DO
1157	END DO
1158
1159	!-----------------------------------------------
1160	! Computation of the snow/ice albedo system
1161	!-------------------------- ---------------------
1162
1163	! Albedo of snow-ice for clear sky.
1164	!-----------------------------------------------
1165	DO jj = 1, jpj
1166	DO ji = 1, jpi
1167	! Case of ice covered by snow.
1168	! melting snow
1169	zihsc1 = 1.0 - MAX ( zzero , SIGN ( zone , - ( hsnif(ji,jj) - c1 ) ) )
1170	zalbpsnm = ( 1.0 - zihsc1 ) &
1171	* ( zficeth(ji,jj) + hsnif(ji,jj) * ( alphd - zficeth(ji,jj) ) / c1 ) &
1172	& + zihsc1 * alphd
1173	! freezing snow
1174	zihsc2 = MAX ( zzero , SIGN ( zone , hsnif(ji,jj) - c2 ) )
1175	zalbpsnf = ( 1.0 - zihsc2 ) * &
1176	( albice + hsnif(ji,jj) * ( alphc - albice ) / c2 ) &
1177	& + zihsc2 * alphc
1178
1179	zitmlsn = MAX ( zzero , SIGN ( zone , sist(ji,jj) - rt0_snow ) )
1180	zalbpsn = zitmlsn * zalbpsnf + ( 1.0 - zitmlsn ) * zalbpsnm
1181
1182	! Case of ice free of snow.
1183	zalbpic = zficeth(ji,jj)
1184
1185	! albedo of the system
1186	zithsn = 1.0 - MAX ( zzero , SIGN ( zone , - hsnif(ji,jj) ) )
1187	palbp(ji,jj) = zithsn * zalbpsn + ( 1.0 - zithsn ) * zalbpic
1188	END DO
1189	END DO
1190
1191	! Albedo of snow-ice for overcast sky.
1192	!----------------------------------------------
1193	palb(:,:) = palbp(:,:) + cgren
1194
1195	!--------------------------------------------
1196	! Computation of the albedo of the ocean
1197	!-------------------------- -----------------
1198
1199
1200	! Parameterization of Briegled and Ramanathan, 1982
1201	zmue14 = zmue**1.4
1202	palcnp(:,:) = 0.05 / ( 1.1 * zmue14 + 0.15 )
1203
1204	! Parameterization of Kondratyev, 1969 and Payne, 1972
1205	palcn(:,:) = 0.06
1206
1207	END SUBROUTINE flx_blk_albedo
1208
1209

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

Download in other formats:

Original Format