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

Context Navigation

source: branches/2013/dev_CMCC_INGV_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90 @ 4217

Last change on this file since 4217 was 4217, checked in by poddo, 11 years ago
Solved problems with namelist parameter land/sea mask
File size: 30.8 KB

Line
1	MODULE sbcblk_mfs
2	!!======================================================================
3	!! * MODULE sbcblk_mfs *
4	!! Ocean forcing: momentum, heat and freshwater flux formulation
5	!!=====================================================================
6	!! History : 3.3 ! 2010-05 (P. Oddo) Original Code
7	!!----------------------------------------------------------------------
8
9	!!----------------------------------------------------------------------
10	!! sbc_blk_mfs : bulk formulation as ocean surface boundary condition
11	!! (forced mode, mfs bulk formulae)
12	!! blk_oce_mfs : ocean: computes momentum, heat and freshwater fluxes
13	!!----------------------------------------------------------------------
14	USE oce ! ocean dynamics and tracers
15	USE dom_oce ! ocean space and time domain
16	USE phycst ! physical constants
17	USE fldread ! read input fields
18	USE sbc_oce ! Surface boundary condition: ocean fields
19	USE iom ! I/O manager library
20	USE in_out_manager ! I/O manager
21	USE lib_mpp ! distribued memory computing library
22	USE wrk_nemo ! work arrays
23	USE timing ! Timing
24	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
25	USE prtctl ! Print control
26	USE sbcwave,ONLY : cdn_wave !wave module
27
28	IMPLICIT NONE
29	PRIVATE
30
31	PUBLIC sbc_blk_mfs ! routine called in sbcmod module
32
33	INTEGER , PARAMETER :: jpfld = 7 ! maximum number of files to read
34	INTEGER , PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point
35	INTEGER , PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point
36	INTEGER , PARAMETER :: jp_clc = 3 ! index of total cloud cover ( % )
37	INTEGER , PARAMETER :: jp_msl = 4 ! index of mean sea level pressure (Pa)
38	INTEGER , PARAMETER :: jp_tair = 5 ! index of 10m air temperature (Kelvin)
39	INTEGER , PARAMETER :: jp_rhm = 6 ! index of dew point temperature (Kelvin)
40	INTEGER , PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s)
41	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
42
43	!! * Substitutions
44	# include "domzgr_substitute.h90"
45	# include "vectopt_loop_substitute.h90"
46	!!----------------------------------------------------------------------
47	!! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
48	!! $Id: sbcblk_mfs.F90 1730 2009-11-16 14:34:19Z poddo $
49	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
50	!!----------------------------------------------------------------------
51
52	CONTAINS
53
54
55	SUBROUTINE sbc_blk_mfs( kt )
56	!!---------------------------------------------------------------------
57	!! * ROUTINE sbc_blk_mfs *
58	!!
59	!! ** Purpose : provide at each time step the surface ocean fluxes
60	!! (momentum, heat, freshwater, runoff is added later in the code)
61	!!
62	!! ** Method : (1) READ Atmospheric data from NetCDF files:
63	!! the 10m wind velocity (i-component) (m/s) at T-point
64	!! the 10m wind velocity (j-component) (m/s) at T-point
65	!! the 2m Dew point Temperature (k)
66	!! the Cloud COver (%)
67	!! the 2m air temperature (Kelvin)
68	!! the Mean Sea Level Preesure (hPa)
69	!! the Climatological Precipitation (kg/m2/s)
70	!! (2) CALL blk_oce_mfs
71	!!
72	!! Computes:
73	!! Solar Radiation using Reed formula (1975, 1977)
74	!! Net Long wave radiation using Bignami et al. (1995)
75	!! Latent and Sensible heat using Kondo (1975)
76	!! Drag coeff using Hllerman and Rosenstein (1983)
77	!! C A U T I O N : never mask the surface stress fields
78	!! the stress is assumed to be in the mesh referential
79	!! i.e. the (i,j) referential
80	!!
81	!! ** Action : defined at each time-step at the air-sea interface
82	!! - utau, vtau i- and j-component of the wind stress
83	!! - taum wind stress module at T-point
84	!! - wndm 10m wind module at T-point
85	!! - qns, qsr non-slor and solar heat flux
86	!! - emp evaporation minus precipitation
87	!!----------------------------------------------------------------------
88	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sh_now ! specific humidity at T-point
89	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: catm ! Cover
90	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: alonl ! Lon
91	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: alatl ! Lat
92	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsst ! SST
93	!!---------------------------------------------------------------------
94	!! Local fluxes variables
95	!!---------------------------------------------------------------------
96	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: qbw ! Net Long wave
97	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ha ! Sesnible
98	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: elat ! Latent
99	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: evap ! evaporation rate
100
101	INTEGER, INTENT( in ) :: kt ! ocean time step
102	!!
103	INTEGER :: ierror ! return error code
104	INTEGER :: ifpr ! dummy loop indice
105	INTEGER :: jj,ji ! dummy loop arguments
106	REAL(wp) :: act_hour
107	!!--------------------------------------------------------------------
108	!! Variables for specific humidity computation
109	!!--------------------------------------------------------------------
110	REAL(wp) :: onsea,par1,par2
111	DATA onsea,par1,par2 / 0.98, 640380., -5107.4 /
112	!! par1 [Kg/m3], par2 [K]
113
114	CHARACTER(len=100) :: cn_dir ! Root directory for location of Atmospheric forcing files
115	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
116	TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_clc, sn_msl ! informations about the fields to be read
117	TYPE(FLD_N) :: sn_tair , sn_rhm, sn_prec ! " "
118	!!---------------------------------------------------------------------
119	NAMELIST/namsbc_mfs/ cn_dir , &
120	& sn_wndi , sn_wndj, sn_clc , sn_msl , &
121	& sn_tair , sn_rhm , sn_prec
122	!!---------------------------------------------------------------------
123	!
124	IF( nn_timing == 1 ) CALL timing_start('sbc_blk_mfs')
125	!
126	! ! ====================== !
127	IF( kt == nit000 ) THEN ! First call kt=nit000 !
128	! ! ====================== !
129	ALLOCATE( sh_now(jpi,jpj), catm(jpi,jpj), alonl(jpi,jpj), alatl(jpi,jpj), &
130	& gsst(jpi,jpj), qbw(jpi,jpj), ha(jpi,jpj), elat(jpi,jpj), &
131	& evap(jpi,jpj), STAT=ierror )
132
133	IF( ierror /= 0 ) CALL ctl_warn('sbc_blk_mfs: failed to allocate arrays')
134	!
135	! set file information (default values)
136	cn_dir = './' ! directory in which the model is executed
137	!
138	! (NB: frequency positive => hours, negative => months)
139	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! land/sea mask !
140	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! filename !
141	sn_wndi = FLD_N( 'ecmwf' , 24 , 'u10' , .false. , .false. , 'daily' , '' , &
142	& '' , '' )
143	sn_wndj = FLD_N( 'ecmwf' , 24 , 'v10' , .false. , .false. , 'daily' , '' , &
144	& '' , '' )
145	sn_clc = FLD_N( 'ecmwf' , 24 , 'clc' , .false. , .false. , 'daily' , '' , &
146	& '' , '' )
147	sn_msl = FLD_N( 'ecmwf' , 24 , 'msl' , .false. , .false. , 'daily' , '' , &
148	& '' , '' )
149	sn_tair = FLD_N( 'ecmwf' , 24 , 't2' , .false. , .false. , 'daily' , '' , &
150	& '' , '' )
151	sn_rhm = FLD_N( 'ecmwf' , 24 , 'rh' , .false. , .false. , 'daily' , '' , &
152	& '' , '' )
153	sn_prec = FLD_N( 'precip_cmap' , -1 , 'precip' , .true. , .true. , 'yearly' , '' , &
154	& '' , '' )
155	!
156	REWIND( numnam ) ! ... read in namlist namsbc_mfs
157	READ ( numnam, namsbc_mfs )
158	!
159	! store namelist information in an array
160	slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj
161	slf_i(jp_clc ) = sn_clc ; slf_i(jp_msl ) = sn_msl
162	slf_i(jp_tair) = sn_tair ; slf_i(jp_rhm) = sn_rhm
163	slf_i(jp_prec) = sn_prec ;
164	!
165	ALLOCATE( sf(jpfld), STAT=ierror ) ! set sf structure
166	IF( ierror > 0 ) THEN
167	CALL ctl_stop( 'sbc_blk_mfs: unable to allocate sf structure' ) ; RETURN
168	ENDIF
169	DO ifpr= 1, jpfld
170	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
171	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
172	END DO
173	! fill sf with slf_i and control print
174	CALL fld_fill( sf, slf_i, cn_dir,'sbc_blk_mfs','bulk formulation for ocean SBC', 'namsbc_mfs' )
175	!
176	ENDIF
177
178	CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step
179
180	catm(:,:) = 0.0 ! initializze cloud cover variable
181	sh_now(:,:) = 0.0 ! initializze specifif humidity variable
182
183	DO jj = 2, jpjm1
184	DO ji = fs_2, fs_jpim1 ! vector opt.
185
186	! Calculate Specific Humidity
187	!-------------------------------------------------
188	sh_now(ji,jj) = (1/1.22) * onsea * par1 * EXP(par2/sf(jp_rhm)%fnow(ji,jj,1))
189
190	! Normalize Clouds
191	!-------------------------------------------------
192	catm(ji,jj) = max(0.0,min(1.0,sf(jp_clc)%fnow(ji,jj,1)*0.01))
193
194	END DO
195	END DO
196
197	! wind module at 10m
198	!--------------------------------------------------
199	wndm(:,:) = SQRT( sf(jp_wndi)%fnow(:,:,1) * sf(jp_wndi)%fnow(:,:,1) &
200	& + sf(jp_wndj)%fnow(:,:,1) * sf(jp_wndj)%fnow(:,:,1) )
201
202	! Some conv for fluxes computation
203	!-------------------------------------------------
204	alonl(:,:) = glamt(:,:) * rad
205	alatl(:,:) = gphit(:,:) * rad
206	gsst(:,:) = tsn(:,:,1,jp_tem) * tmask(:,:,1)
207
208	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
209
210	! Force to zero the output of fluxes
211	!-------------------------------------------------
212	qsr(:,:) = 0.0 ; qbw(:,:) = 0.0 ;
213	ha(:,:) = 0.0 ; elat(:,:) = 0.0 ;
214	evap(:,:) = 0.0 ; utau(:,:) = 0.0 ;
215	vtau(:,:) = 0.0
216
217	CALL lbc_lnk( sf(jp_wndi)%fnow(:,:,1), 'T', -1. )
218	CALL lbc_lnk( sf(jp_wndj)%fnow(:,:,1), 'T', -1. )
219
220	act_hour = (( nsec_year / rday ) - INT (nsec_year / rday)) * rjjhh
221
222	CALL fluxes_mfs(alatl,alonl,act_hour, & ! input static
223	gsst(:,:),sf(jp_tair)%fnow(:,:,1),sh_now(:,:), & ! input dynamic
224	sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), & ! input dynamic
225	sf(jp_msl)%fnow(:,:,1) , catm(:,:) , & ! input dynamic
226	qsr,qbw,ha,elat,evap,utau,vtau) ! output
227
228	! Shortwave radiation
229	!--------------------------------------------------
230	qsr(:,:) = qsr(:,:) * tmask(:,:,1)
231
232	! total non solar heat flux over water
233	!--------------------------------------------------
234	qns(:,:) = -1. * ( qbw(:,:) + ha(:,:) + elat(:,:) )
235	qns(:,:) = qns(:,:)*tmask(:,:,1)
236
237	! mask the wind module at 10m
238	!--------------------------------------------------
239	wndm(:,:) = wndm(:,:) * tmask(:,:,1)
240
241	! wind stress module (taum) into T-grid
242	!--------------------------------------------------
243	taum(:,:) = SQRT( utau(:,:) * utau(:,:) + vtau(:,:) * vtau(:,:) ) * tmask(:,:,1)
244
245	CALL lbc_lnk( taum, 'T', 1. )
246
247	! Interpolate utau, vtau into the grid_V and grid_V
248	!-------------------------------------------------
249
250	DO jj = 1, jpjm1
251	DO ji = 1, fs_jpim1
252	utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( utau(ji,jj) * tmask(ji,jj,1) &
253	& + utau(ji+1,jj) * tmask(ji+1,jj,1) )
254	vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( vtau(ji,jj) * tmask(ji,jj,1) &
255	& + vtau(ji,jj+1) * tmask(ji,jj+1,1) )
256	END DO
257	END DO
258
259	CALL lbc_lnk( utau(:,:), 'U', -1. )
260	CALL lbc_lnk( vtau(:,:), 'V', -1. )
261
262	! for basin budget and cooerence
263	!--------------------------------------------------
264	!CDIR COLLAPSE
265	emp (:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1)
266	!CDIR COLLAPSE
267
268	CALL iom_put( "qlw_oce", qbw ) ! output downward longwave heat over the ocean
269	CALL iom_put( "qsb_oce", - ha ) ! output downward sensible heat over the ocean
270	CALL iom_put( "qla_oce", - elat ) ! output downward latent heat over the ocean
271	CALL iom_put( "qns_oce", qns ) ! output downward non solar heat over the ocean
272
273	ENDIF
274	!
275	IF( nn_timing == 1 ) CALL timing_stop('sbc_blk_mfs')
276	!
277	END SUBROUTINE sbc_blk_mfs
278
279
280	SUBROUTINE fluxes_mfs(alat,alon,hour, &
281	sst,tnow,shnow,unow,vnow,mslnow,cldnow,qsw,qbw,ha,elat, &
282	evap,taux,tauy)
283	!!----------------------------------------------------------------------
284	!! * ROUTINE fluxes_mfs *
285	!!
286	!! --- it provides SURFACE HEAT and MOMENTUM FLUXES in MKS :
287	!!
288	!! 1) Water flux (WFLUX) [ watt/m*m ]
289	!! 2) Short wave flux (QSW) [ watt/m*m ] Reed 1977
290	!! 3) Long wave flux backward (QBW) [ watt/m*m ]
291	!! 4) Latent heat of evaporation (ELAT) [ watt/m*m ]
292	!! 5) Sensible heat flux (HA) [ watt/m*m ]
293	!! 6) Wind stress x-component (TAUX) [ newton/m*m ]
294	!! 7) Wind stress y-component (TAUY) [ newton/m*m ]
295	!!
296	!!----------------------------------------------------------------------
297	USE sbcblk_core, ONLY: turb_core_2z ! For wave coupling and Tair/rh from 2 to 10m
298
299	REAL(wp), INTENT(in ) :: hour
300	REAL(wp), INTENT(in ), DIMENSION (:,:) :: sst, unow, alat , alon
301	REAL(wp), INTENT(in ), DIMENSION (:,:) :: vnow, cldnow, mslnow
302	REAL(wp), INTENT(out ), DIMENSION (:,:) :: qsw, qbw, ha, elat
303	REAL(wp), INTENT(out ), DIMENSION (:,:) :: evap,taux,tauy
304	REAL(wp), INTENT(inout), DIMENSION (:,:) :: tnow , shnow
305
306	INTEGER :: ji,jj
307	REAL(wp) :: wair, vtnow, ea, deltemp, s, stp , fh , fe
308	REAL(wp) :: esre, cseep
309
310	REAL(wp), DIMENSION (:,:), POINTER :: rspeed, sh10now, t10now, cdx, ce, shms
311	REAL(wp), DIMENSION (:,:), POINTER :: rhom, sstk, ch, rel_windu, rel_windv
312	!!----------------------------------------------------------------------
313	!! coefficients ( in MKS ) :
314	!!----------------------------------------------------------------------
315
316	REAL(wp), PARAMETER :: ps = 1013. ! --- surface air pressure
317	REAL(wp), PARAMETER :: expsi=0.622 ! --- expsi
318	REAL(wp), PARAMETER :: rd=287. ! --- dry air gas constant
319	REAL(wp), PARAMETER :: cp=1005. ! --- specific heat capacity
320	REAL(wp), PARAMETER :: onsea=0.98 ! --- specific humidity factors
321	REAL(wp), PARAMETER :: par1=640380. ! [Kg/m3]
322	REAL(wp), PARAMETER :: par2=-5107.4 ! [K]
323
324	!---------------------------------------------------------------------
325	!--- define Kondo parameters
326	!---------------------------------------------------------------------
327
328	REAL(wp), DIMENSION(5) :: a_h = (/0.0,0.927,1.15,1.17,1.652/)
329	REAL(wp), DIMENSION(5) :: a_e = (/0.0,0.969,1.18,1.196,1.68/)
330	REAL(wp), DIMENSION(5) :: b_h = (/1.185,0.0546,0.01,0.0075,-0.017/)
331	REAL(wp), DIMENSION(5) :: b_e = (/1.23,0.0521,0.01,0.008,-0.016/)
332	REAL(wp), DIMENSION(5) :: c_h = (/0.0,0.0,0.0,-0.00045,0.0/)
333	REAL(wp), DIMENSION(5) :: c_e = (/0.0,0.0,0.0,-0.0004,0.0/)
334	REAL(wp), DIMENSION(5) :: p_h = (/-0.157,1.0,1.0,1.0,1.0/)
335	REAL(wp), DIMENSION(5) :: p_e = (/-0.16,1.0,1.0,1.0,1.0/)
336	INTEGER :: kku !index varing with wind speed
337	!
338	IF( nn_timing == 1 ) CALL timing_start('fluxes_mfs')
339	!
340	CALL wrk_alloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
341	CALL wrk_alloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )
342
343	!!----------------------------------------------------------------------
344	!! ------------------ (i) short wave
345	!!----------------------------------------------------------------------
346
347	CALL qshort(hour,alat,alon,cldnow,qsw)
348
349	rel_windu(:,:) = 0.0_wp
350	rel_windv(:,:) = 0.0_wp
351
352	DO jj = 2, jpj
353	DO ji = fs_2, jpi
354	rel_windu(ji,jj) = unow(ji,jj) - 0.5_wp * ( ssu_m(ji-1,jj) + ssu_m(ji,jj) )
355	rel_windv(ji,jj) = vnow(ji,jj) - 0.5_wp * ( ssv_m(ji,jj-1) + ssv_m(ji,jj) )
356	END DO
357	END DO
358
359	rspeed(:,:)= SQRT(rel_windu(:,:)*rel_windu(:,:) &
360	& + rel_windv(:,:)*rel_windv(:,:))
361
362	sstk(:,:) = sst(:,:) + rtt !- SST data in Kelvin degrees
363	shms(:,:) = (1/1.22)onseapar1*EXP(par2/sstk(:,:)) !- Saturation Specific Humidity
364
365	! --- Transport temperature and humidity from 2m to 10m
366	!----------------------------------------------------------------------
367
368	t10now(:,:) = 0.0 ; sh10now(:,:)= 0.0
369	! Note that air temp is converted in potential temp
370	CALL turb_core_2z(2.,10.,sstk,tnow+2*0.0098,shms,shnow,rspeed, &
371	& Cdx,Ch,Ce,t10now,sh10now )
372	tnow(:,:) = t10now(:,:) ; shnow(:,:) = sh10now(:,:)
373
374	!!----------------------------------------------------------------------
375	!! ------------------ (ii) net long wave
376	!!----------------------------------------------------------------------
377
378	DO jj = 1, jpj
379	DO ji = 1, jpi
380	wair = shnow(ji,jj) / (1 - shnow(ji,jj)) ! mixing ratio of the air (Wallace and Hobbs)
381	vtnow = (tnow(ji,jj)(expsi+wair))/(expsi(1.+wair)) ! virtual temperature of air
382	rhom(ji,jj) = 100.*(ps/rd)/vtnow ! density of the moist air
383	ea = (wair / (wair + 0.622 )) * mslnow(ji,jj)
384
385	qbw(ji,jj) = emicstefan( sstk(ji,jj)**4. ) &
386	- ( stefan( tnow(ji,jj)4. ) ( 0.653 + 0.00535*ea ) ) &
387	* ( 1. + 0.1762( cldnow(ji,jj)*2. ) )
388
389	END DO
390	END DO
391
392	DO jj = 1, jpj
393	DO ji = 1, jpi
394	!!----------------------------------------------------------------------
395	!! ------------------ (iii) sensible heat
396	!!----------------------------------------------------------------------
397
398	!! --- calculates the term : ( Ts - Ta )
399	!!----------------------------------------------------------------------
400	deltemp = sstk(ji,jj) - tnow (ji,jj)
401
402	!!----------------------------------------------------------------------
403	!! --- variable turbulent exchange coefficients ( from Kondo 1975 )
404	!! --- calculate the Neutral Transfer Coefficent using an empiric formula
405	!! --- by Kondo et al. Then it applies the diabatic approximation.
406	!!----------------------------------------------------------------------
407
408	s = deltemp/(wndm(ji,jj)**2.) !! --- calculate S
409	stp = s*abs(s)/(abs(s)+0.01) !! --- calculate the Stability Parameter
410
411	!!----------------------------------------------------------------------
412	!! --- for stable condition (sst-t_air < 0):
413	!!----------------------------------------------------------------------
414
415	IF (s.lt.0. .and. ((stp.gt.-3.3).and.(stp.lt.0.))) THEN
416	fh = 0.1_wp+0.03_wpstp+0.9_wpexp(4.8_wp*stp)
417	fe = fh
418	ELSE IF (s.lt.0. .and. stp.le.-3.3) THEN
419	fh = 0._wp
420	fe = fh
421	ELSE ! --- for unstable condition
422	fh = 1.0_wp+0.63_wp*sqrt(stp)
423	fe = fh
424	ENDIF
425
426	!!----------------------------------------------------------------------
427	!! --- calculate the coefficient CH,CE,CD
428	!!----------------------------------------------------------------------
429
430	IF (wndm(ji,jj) >= 0. .AND. wndm(ji,jj) <= 2.2) THEN
431	kku=1
432	ELSE IF (wndm(ji,jj) > 2.2 .AND. wndm(ji,jj) <= 5.0) THEN
433	kku=2
434	ELSE IF (wndm(ji,jj) > 5.0 .AND. wndm(ji,jj) <= 8.0) THEN
435	kku=3
436	ELSE IF (wndm(ji,jj) > 8.0 .AND. wndm(ji,jj) <= 25.0) THEN
437	kku=4
438	ELSE IF (wndm(ji,jj) > 25.0 ) THEN
439	kku=5
440	ENDIF
441
442	ch(ji,jj) = ( a_h(kku) + b_h(kku) * wndm(ji,jj) ** p_h(kku) &
443	+ c_h(kku) * (wndm(ji,jj)-8 ) *2) fh
444
445	ce(ji,jj) = ( a_e(kku) + b_e(kku) * wndm(ji,jj) ** p_e(kku) &
446	+ c_e(kku) * (wndm(ji,jj)-8 ) *2) fe
447
448	ch(ji,jj) = ch(ji,jj) / 1000.0
449	ce(ji,jj) = ce(ji,jj) / 1000.0
450
451	IF (wndm(ji,jj)<0.3) THEN
452	ch(ji,jj) = 1.3e-03 * fh
453	ce(ji,jj) = 1.5e-03 * fe
454	ELSE IF(wndm(ji,jj)>50.0) THEN
455	ch(ji,jj) = 1.25e-03 * fh
456	ce(ji,jj) = 1.30e-03 * fe
457	ENDIF
458
459	!!----------------------------------------------------------------------
460	!! --- calculates the SENSIBLE HEAT FLUX in MKS ( watt/m*m )
461	!!----------------------------------------------------------------------
462
463	HA(ji,jj) = rhom(ji,jj)cpch(ji,jj)wndm(ji,jj)deltemp
464
465	!!----------------------------------------------------------------------
466	!! ------------------ (iv) latent heat
467	!! --- calculates the LATENT HEAT FLUX ( watt/m*m )
468	!! --- ELAT = LrhoCe\|V\|[qs(Ts)-qa(t2d)]
469	!!----------------------------------------------------------------------
470
471	esre = shms(ji,jj) - shnow(ji,jj) ! --- calculates the term : qs(Ta)-qa(t2d)
472
473	cseep = ce(ji,jj) * wndm(ji,jj) * esre ! --- calculates the term : Ce\|V\|[qs(Ts)-qa(t2d)]
474
475	evap(ji,jj) = (cseep * rhom(ji,jj)) ! in [kg/m2/sec] !! --- calculates the EVAPORATION RATE [m/yr]
476
477	elat(ji,jj) = rhom(ji,jj) * cseep * heatlat(sst(ji,jj))
478
479	!!----------------------------------------------------------------------
480	!! --- calculates the Drag Coefficient
481	!!----------------------------------------------------------------------
482
483	!!----------------------------------------------------------------------
484	!! --- deltemp should be (Ts - Ta) in the formula estimating
485	!! --- drag coefficient
486	!!----------------------------------------------------------------------
487
488	IF( .NOT. ln_cdgw ) THEN
489	cdx(ji,jj) = cd_HR(wndm(ji,jj),deltemp)
490	ENDIF
491
492	END DO
493	END DO
494
495	!!----------------------------------------------------------------------
496	!! --- calculates the wind stresses in MKS ( newton/m*m )
497	!! --- taux= rhoCd\|V\|u tauy= rhoCd\|V\|v
498	!!----------------------------------------------------------------------
499
500	taux(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windu(:,:)
501	tauy(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windv(:,:)
502
503	CALL wrk_dealloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
504	CALL wrk_dealloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )
505	!
506	IF( nn_timing == 1 ) CALL timing_stop('fluxes_mfs')
507	!
508	END SUBROUTINE fluxes_mfs
509
510
511	REAL(wp) FUNCTION cd_HR(speed,delt)
512	!!----------------------------------------------------------------------
513	!! --- calculates the Drag Coefficient as a function of the abs. value
514	!! --- of the wind velocity ( Hellermann and Rosenstein )
515	!!----------------------------------------------------------------------
516
517	REAL(wp), INTENT(in) :: speed,delt
518	REAL(wp), PARAMETER :: a1=0.934e-3 , a2=0.788e-4, a3=0.868e-4
519	REAL(wp), PARAMETER :: a4=-0.616e-6, a5=-.120e-5, a6=-.214e-5
520
521	cd_HR = a1 + a2speed + a3delt + a4speedspeed &
522	+ a5deltdelt + a6speeddelt
523
524	END FUNCTION cd_HR
525
526	REAL(wp) function HEATLAT(t)
527	!!----------------------------------------------------------------------
528	!! --- calculates the Latent Heat of Vaporization ( J/kg ) as function of
529	!! --- the temperature ( Celsius degrees )
530	!! --- ( from A. Gill pag. 607 )
531	!!
532	!! --- Constant Latent Heat of Vaporization ( Rosati,Miyakoda 1988 )
533	!! L = 2.501e+6 (MKS)
534	!!----------------------------------------------------------------------
535
536	REAL(wp) , intent(in) :: t
537
538	heatlat = 2.5008e+6 -2.3e+3*t
539
540	END FUNCTION HEATLAT
541
542
543	SUBROUTINE qshort(hour,alat,alon,cldnow,qsw)
544	!!----------------------------------------------------------------------
545	!! * ROUTINE qshort *
546	!!
547	!! ** Purpose : Compute Solar Radiation
548	!!
549	!! ** Method : Compute Solar Radiation according Astronomical
550	!! formulae
551	!!
552	!! References : Reed RK (1975) and Reed RK (1977)
553	!!
554	!! Note: alat,alon - (lat, lon) in radians
555	!!----------------------------------------------------------------------
556	REAL(wp), INTENT (in) :: hour
557
558	REAL(wp), INTENT(in ), DIMENSION(:,:) :: alat,alon
559	REAL(wp), INTENT(in ), DIMENSION(:,:) :: cldnow
560	REAL(wp), INTENT(out), DIMENSION(:,:) :: qsw
561	REAL(wp), DIMENSION(12) :: alpham
562
563	REAL(wp), PARAMETER :: eclips=23.439* (3.141592653589793_wp / 180._wp)
564	REAL(wp), PARAMETER :: solar = 1350.
565	REAL(wp), PARAMETER :: tau = 0.7
566	REAL(wp), PARAMETER :: aozone = 0.09
567	REAL(wp), PARAMETER :: yrdays = 360.
568	REAL(wp) :: days, th0,th02,th03, sundec, thsun, coszen, qatten
569	REAL(wp) :: qzer, qdir,qdiff,qtot,tjul,sunbet
570	REAL(wp) :: albedo
571	INTEGER :: jj, ji
572
573	!!----------------------------------------------------------------------
574	!! --- albedo monthly values from Payne (1972) as means of the values
575	!! --- at 40N and 30N for the Atlantic Ocean ( hence the same latitudinal
576	!! --- band of the Mediterranean Sea ) :
577	!!----------------------------------------------------------------------
578
579	data alpham /0.095,0.08,0.065,0.065,0.06,0.06,0.06,0.06, &
580	0.065,0.075,0.09,0.10/
581
582	!!----------------------------------------------------------------------
583	!! days is the number of days elapsed until the day=nday_year
584	!!----------------------------------------------------------------------
585	days = nday_year -1.
586	th0 = 2.rpidays/yrdays
587	th02 = 2.*th0
588	th03 = 3.*th0
589
590	!! --- sun declination :
591	!!----------------------------------------------------------------------
592	sundec = 0.006918 - 0.399912cos(th0) + 0.070257sin(th0) - &
593	0.006758cos(th02) + 0.000907sin(th02) - &
594	0.002697cos(th03) + 0.001480sin(th03)
595
596	DO jj = 1, jpj
597	DO ji = 1, jpi
598
599	!! --- sun hour angle :
600	!!----------------------------------------------------------------------
601	thsun = (hour -12.)15.rad + alon(ji,jj)
602
603	!! --- cosine of the solar zenith angle :
604	!!----------------------------------------------------------------------
605	coszen =sin(alat(ji,jj))*sin(sundec) &
606	+cos(alat(ji,jj))cos(sundec)cos(thsun)
607
608	IF(coszen .LE. 5.035D-04) THEN
609	coszen = 0.0
610	qatten = 0.0
611	ELSE
612	qatten = tau**(1./coszen)
613	END IF
614
615	qzer = coszen * solar * &
616	(1.0+1.67E-2cos(rpi2.(days-3.0)/365.0))*2
617	qdir = qzer * qatten
618	qdiff = ((1.-aozone)qzer - qdir) 0.5
619	qtot = qdir + qdiff
620	tjul = (days -81.)*rad
621
622	!! --- sin of the solar noon altitude in radians :
623	!!----------------------------------------------------------------------
624	sunbet=sin(alat(ji,jj))sin(eclipssin(tjul)) + &
625	cos(alat(ji,jj))cos(eclipssin(tjul))
626
627	!! --- solar noon altitude in degrees :
628	!!----------------------------------------------------------------------
629
630	sunbet = asin(sunbet)/rad
631
632	!!----------------------------------------------------------------------
633	!! --- calculates the albedo according to Payne (1972)
634	!!----------------------------------------------------------------------
635
636	albedo = alpham(nmonth)
637
638	!!----------------------------------------------------------------------
639	!! --- ( radiation as from Reed(1977), Simpson and Paulson(1979) )
640	!! --- calculates SHORT WAVE FLUX ( watt/m*m )
641	!! --- ( Rosati,Miyakoda 1988 ; eq. 3.8 )
642	!!----------------------------------------------------------------------
643
644	IF(cldnow(ji,jj).LT.0.3) THEN
645	qsw(ji,jj) = qtot * (1.-albedo)
646	ELSE
647	qsw(ji,jj) = qtot(1.-0.62cldnow(ji,jj) &
648	+ .0019sunbet)(1.-albedo)
649	ENDIF
650
651	END DO
652	END DO
653
654	END SUBROUTINE qshort
655
656
657	!!======================================================================
658
659	END MODULE sbcblk_mfs

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

Download in other formats: