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

Context Navigation

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90 @ 3152

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

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

Download in other formats: