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 @ 473

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

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

Download in other formats:

Original Format