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/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90 @ 3402

Last change on this file since 3402 was 3402, checked in by acc, 12 years ago
Branch: dev_r3385_NOCS04_HAMF; #665. Stage 2 of 2012 development: suppression of emps array and introduction of sfx (salt flux) array with associated code to setup the options for embedding the seaice into the ocean
File size: 30.5 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 !
140	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs !
141	sn_wndi = FLD_N( 'ecmwf' , 24 , 'u10' , .false. , .false. , 'daily' , '' , '' )
142	sn_wndj = FLD_N( 'ecmwf' , 24 , 'v10' , .false. , .false. , 'daily' , '' , '' )
143	sn_clc = FLD_N( 'ecmwf' , 24 , 'clc' , .false. , .false. , 'daily' , '' , '' )
144	sn_msl = FLD_N( 'ecmwf' , 24 , 'msl' , .false. , .false. , 'daily' , '' , '' )
145	sn_tair = FLD_N( 'ecmwf' , 24 , 't2' , .false. , .false. , 'daily' , '' , '' )
146	sn_rhm = FLD_N( 'ecmwf' , 24 , 'rh' , .false. , .false. , 'daily' , '' , '' )
147	sn_prec = FLD_N( 'precip_cmap' , -1 , 'precip' , .true. , .true. , 'yearly' , '' , '' )
148	!
149	REWIND( numnam ) ! ... read in namlist namsbc_mfs
150	READ ( numnam, namsbc_mfs )
151	!
152	! store namelist information in an array
153	slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj
154	slf_i(jp_clc ) = sn_clc ; slf_i(jp_msl ) = sn_msl
155	slf_i(jp_tair) = sn_tair ; slf_i(jp_rhm) = sn_rhm
156	slf_i(jp_prec) = sn_prec ;
157	!
158	ALLOCATE( sf(jpfld), STAT=ierror ) ! set sf structure
159	IF( ierror > 0 ) THEN
160	CALL ctl_stop( 'sbc_blk_mfs: unable to allocate sf structure' ) ; RETURN
161	ENDIF
162	DO ifpr= 1, jpfld
163	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
164	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
165	END DO
166	! fill sf with slf_i and control print
167	CALL fld_fill( sf, slf_i, cn_dir,'sbc_blk_mfs','bulk formulation for ocean SBC', 'namsbc_mfs' )
168	!
169	ENDIF
170
171	CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step
172
173	catm(:,:) = 0.0 ! initializze cloud cover variable
174	sh_now(:,:) = 0.0 ! initializze specifif humidity variable
175
176	DO jj = 2, jpjm1
177	DO ji = fs_2, fs_jpim1 ! vector opt.
178
179	! Calculate Specific Humidity
180	!-------------------------------------------------
181	sh_now(ji,jj) = (1/1.22) * onsea * par1 * EXP(par2/sf(jp_rhm)%fnow(ji,jj,1))
182
183	! Normalize Clouds
184	!-------------------------------------------------
185	catm(ji,jj) = max(0.0,min(1.0,sf(jp_clc)%fnow(ji,jj,1)*0.01))
186
187	END DO
188	END DO
189
190	! wind module at 10m
191	!--------------------------------------------------
192	wndm(:,:) = SQRT( sf(jp_wndi)%fnow(:,:,1) * sf(jp_wndi)%fnow(:,:,1) &
193	& + sf(jp_wndj)%fnow(:,:,1) * sf(jp_wndj)%fnow(:,:,1) )
194
195	! Some conv for fluxes computation
196	!-------------------------------------------------
197	alonl(:,:) = glamt(:,:) * rad
198	alatl(:,:) = gphit(:,:) * rad
199	gsst(:,:) = tsn(:,:,1,jp_tem) * tmask(:,:,1)
200
201	IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
202
203	! Force to zero the output of fluxes
204	!-------------------------------------------------
205	qsr(:,:) = 0.0 ; qbw(:,:) = 0.0 ;
206	ha(:,:) = 0.0 ; elat(:,:) = 0.0 ;
207	evap(:,:) = 0.0 ; utau(:,:) = 0.0 ;
208	vtau(:,:) = 0.0
209
210	CALL lbc_lnk( sf(jp_wndi)%fnow(:,:,1), 'T', -1. )
211	CALL lbc_lnk( sf(jp_wndj)%fnow(:,:,1), 'T', -1. )
212
213	act_hour = (( nsec_year / rday ) - INT (nsec_year / rday)) * rjjhh
214
215	CALL fluxes_mfs(alatl,alonl,act_hour, & ! input static
216	gsst(:,:),sf(jp_tair)%fnow(:,:,1),sh_now(:,:), & ! input dynamic
217	sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), & ! input dynamic
218	sf(jp_msl)%fnow(:,:,1) , catm(:,:) , & ! input dynamic
219	qsr,qbw,ha,elat,evap,utau,vtau) ! output
220
221	! Shortwave radiation
222	!--------------------------------------------------
223	qsr(:,:) = qsr(:,:) * tmask(:,:,1)
224
225	! total non solar heat flux over water
226	!--------------------------------------------------
227	qns(:,:) = -1. * ( qbw(:,:) + ha(:,:) + elat(:,:) )
228	qns(:,:) = qns(:,:)*tmask(:,:,1)
229
230	! mask the wind module at 10m
231	!--------------------------------------------------
232	wndm(:,:) = wndm(:,:) * tmask(:,:,1)
233
234	! wind stress module (taum) into T-grid
235	!--------------------------------------------------
236	taum(:,:) = SQRT( utau(:,:) * utau(:,:) + vtau(:,:) * vtau(:,:) ) * tmask(:,:,1)
237
238	CALL lbc_lnk( taum, 'T', 1. )
239
240	! Interpolate utau, vtau into the grid_V and grid_V
241	!-------------------------------------------------
242
243	DO jj = 1, jpjm1
244	DO ji = 1, fs_jpim1
245	utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( utau(ji,jj) * tmask(ji,jj,1) &
246	& + utau(ji+1,jj) * tmask(ji+1,jj,1) )
247	vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( vtau(ji,jj) * tmask(ji,jj,1) &
248	& + vtau(ji,jj+1) * tmask(ji,jj+1,1) )
249	END DO
250	END DO
251
252	CALL lbc_lnk( utau(:,:), 'U', -1. )
253	CALL lbc_lnk( vtau(:,:), 'V', -1. )
254
255	! for basin budget and cooerence
256	!--------------------------------------------------
257	!CDIR COLLAPSE
258	emp (:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1)
259	!CDIR COLLAPSE
260
261	CALL iom_put( "qlw_oce", qbw ) ! output downward longwave heat over the ocean
262	CALL iom_put( "qsb_oce", - ha ) ! output downward sensible heat over the ocean
263	CALL iom_put( "qla_oce", - elat ) ! output downward latent heat over the ocean
264	CALL iom_put( "qns_oce", qns ) ! output downward non solar heat over the ocean
265
266	ENDIF
267	!
268	IF( nn_timing == 1 ) CALL timing_stop('sbc_blk_mfs')
269	!
270	END SUBROUTINE sbc_blk_mfs
271
272
273	SUBROUTINE fluxes_mfs(alat,alon,hour, &
274	sst,tnow,shnow,unow,vnow,mslnow,cldnow,qsw,qbw,ha,elat, &
275	evap,taux,tauy)
276	!!----------------------------------------------------------------------
277	!! * ROUTINE fluxes_mfs *
278	!!
279	!! --- it provides SURFACE HEAT and MOMENTUM FLUXES in MKS :
280	!!
281	!! 1) Water flux (WFLUX) [ watt/m*m ]
282	!! 2) Short wave flux (QSW) [ watt/m*m ] Reed 1977
283	!! 3) Long wave flux backward (QBW) [ watt/m*m ]
284	!! 4) Latent heat of evaporation (ELAT) [ watt/m*m ]
285	!! 5) Sensible heat flux (HA) [ watt/m*m ]
286	!! 6) Wind stress x-component (TAUX) [ newton/m*m ]
287	!! 7) Wind stress y-component (TAUY) [ newton/m*m ]
288	!!
289	!!----------------------------------------------------------------------
290	USE sbcblk_core, ONLY: turb_core_2z ! For wave coupling and Tair/rh from 2 to 10m
291
292	REAL(wp), INTENT(in ) :: hour
293	REAL(wp), INTENT(in ), DIMENSION (:,:) :: sst, unow, alat , alon
294	REAL(wp), INTENT(in ), DIMENSION (:,:) :: vnow, cldnow, mslnow
295	REAL(wp), INTENT(out ), DIMENSION (:,:) :: qsw, qbw, ha, elat
296	REAL(wp), INTENT(out ), DIMENSION (:,:) :: evap,taux,tauy
297	REAL(wp), INTENT(inout), DIMENSION (:,:) :: tnow , shnow
298
299	INTEGER :: ji,jj
300	REAL(wp) :: wair, vtnow, ea, deltemp, s, stp , fh , fe
301	REAL(wp) :: esre, cseep
302
303	REAL(wp), DIMENSION (:,:), POINTER :: rspeed, sh10now, t10now, cdx, ce, shms
304	REAL(wp), DIMENSION (:,:), POINTER :: rhom, sstk, ch, rel_windu, rel_windv
305	!!----------------------------------------------------------------------
306	!! coefficients ( in MKS ) :
307	!!----------------------------------------------------------------------
308
309	REAL(wp), PARAMETER :: ps = 1013. ! --- surface air pressure
310	REAL(wp), PARAMETER :: expsi=0.622 ! --- expsi
311	REAL(wp), PARAMETER :: rd=287. ! --- dry air gas constant
312	REAL(wp), PARAMETER :: cp=1005. ! --- specific heat capacity
313	REAL(wp), PARAMETER :: onsea=0.98 ! --- specific humidity factors
314	REAL(wp), PARAMETER :: par1=640380. ! [Kg/m3]
315	REAL(wp), PARAMETER :: par2=-5107.4 ! [K]
316
317	!---------------------------------------------------------------------
318	!--- define Kondo parameters
319	!---------------------------------------------------------------------
320
321	REAL(wp), DIMENSION(5) :: a_h = (/0.0,0.927,1.15,1.17,1.652/)
322	REAL(wp), DIMENSION(5) :: a_e = (/0.0,0.969,1.18,1.196,1.68/)
323	REAL(wp), DIMENSION(5) :: b_h = (/1.185,0.0546,0.01,0.0075,-0.017/)
324	REAL(wp), DIMENSION(5) :: b_e = (/1.23,0.0521,0.01,0.008,-0.016/)
325	REAL(wp), DIMENSION(5) :: c_h = (/0.0,0.0,0.0,-0.00045,0.0/)
326	REAL(wp), DIMENSION(5) :: c_e = (/0.0,0.0,0.0,-0.0004,0.0/)
327	REAL(wp), DIMENSION(5) :: p_h = (/-0.157,1.0,1.0,1.0,1.0/)
328	REAL(wp), DIMENSION(5) :: p_e = (/-0.16,1.0,1.0,1.0,1.0/)
329	INTEGER :: kku !index varing with wind speed
330	!
331	IF( nn_timing == 1 ) CALL timing_start('fluxes_mfs')
332	!
333	CALL wrk_alloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
334	CALL wrk_alloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )
335
336	!!----------------------------------------------------------------------
337	!! ------------------ (i) short wave
338	!!----------------------------------------------------------------------
339
340	CALL qshort(hour,alat,alon,cldnow,qsw)
341
342	rel_windu(:,:) = 0.0_wp
343	rel_windv(:,:) = 0.0_wp
344
345	DO jj = 2, jpj
346	DO ji = fs_2, jpi
347	rel_windu(ji,jj) = unow(ji,jj) - 0.5_wp * ( ssu_m(ji-1,jj) + ssu_m(ji,jj) )
348	rel_windv(ji,jj) = vnow(ji,jj) - 0.5_wp * ( ssv_m(ji,jj-1) + ssv_m(ji,jj) )
349	END DO
350	END DO
351
352	rspeed(:,:)= SQRT(rel_windu(:,:)*rel_windu(:,:) &
353	& + rel_windv(:,:)*rel_windv(:,:))
354
355	sstk(:,:) = sst(:,:) + rtt !- SST data in Kelvin degrees
356	shms(:,:) = (1/1.22)onseapar1*EXP(par2/sstk(:,:)) !- Saturation Specific Humidity
357
358	! --- Transport temperature and humidity from 2m to 10m
359	!----------------------------------------------------------------------
360
361	t10now(:,:) = 0.0 ; sh10now(:,:)= 0.0
362	! Note that air temp is converted in potential temp
363	CALL turb_core_2z(2.,10.,sstk,tnow+2*0.0098,shms,shnow,rspeed, &
364	& Cdx,Ch,Ce,t10now,sh10now )
365	tnow(:,:) = t10now(:,:) ; shnow(:,:) = sh10now(:,:)
366
367	!!----------------------------------------------------------------------
368	!! ------------------ (ii) net long wave
369	!!----------------------------------------------------------------------
370
371	DO jj = 1, jpj
372	DO ji = 1, jpi
373	wair = shnow(ji,jj) / (1 - shnow(ji,jj)) ! mixing ratio of the air (Wallace and Hobbs)
374	vtnow = (tnow(ji,jj)(expsi+wair))/(expsi(1.+wair)) ! virtual temperature of air
375	rhom(ji,jj) = 100.*(ps/rd)/vtnow ! density of the moist air
376	ea = (wair / (wair + 0.622 )) * mslnow(ji,jj)
377
378	qbw(ji,jj) = emicstefan( sstk(ji,jj)**4. ) &
379	- ( stefan( tnow(ji,jj)4. ) ( 0.653 + 0.00535*ea ) ) &
380	* ( 1. + 0.1762( cldnow(ji,jj)*2. ) )
381
382	END DO
383	END DO
384
385	DO jj = 1, jpj
386	DO ji = 1, jpi
387	!!----------------------------------------------------------------------
388	!! ------------------ (iii) sensible heat
389	!!----------------------------------------------------------------------
390
391	!! --- calculates the term : ( Ts - Ta )
392	!!----------------------------------------------------------------------
393	deltemp = sstk(ji,jj) - tnow (ji,jj)
394
395	!!----------------------------------------------------------------------
396	!! --- variable turbulent exchange coefficients ( from Kondo 1975 )
397	!! --- calculate the Neutral Transfer Coefficent using an empiric formula
398	!! --- by Kondo et al. Then it applies the diabatic approximation.
399	!!----------------------------------------------------------------------
400
401	s = deltemp/(wndm(ji,jj)**2.) !! --- calculate S
402	stp = s*abs(s)/(abs(s)+0.01) !! --- calculate the Stability Parameter
403
404	!!----------------------------------------------------------------------
405	!! --- for stable condition (sst-t_air < 0):
406	!!----------------------------------------------------------------------
407
408	IF (s.lt.0. .and. ((stp.gt.-3.3).and.(stp.lt.0.))) THEN
409	fh = 0.1_wp+0.03_wpstp+0.9_wpexp(4.8_wp*stp)
410	fe = fh
411	ELSE IF (s.lt.0. .and. stp.le.-3.3) THEN
412	fh = 0._wp
413	fe = fh
414	ELSE ! --- for unstable condition
415	fh = 1.0_wp+0.63_wp*sqrt(stp)
416	fe = fh
417	ENDIF
418
419	!!----------------------------------------------------------------------
420	!! --- calculate the coefficient CH,CE,CD
421	!!----------------------------------------------------------------------
422
423	IF (wndm(ji,jj) >= 0. .AND. wndm(ji,jj) <= 2.2) THEN
424	kku=1
425	ELSE IF (wndm(ji,jj) > 2.2 .AND. wndm(ji,jj) <= 5.0) THEN
426	kku=2
427	ELSE IF (wndm(ji,jj) > 5.0 .AND. wndm(ji,jj) <= 8.0) THEN
428	kku=3
429	ELSE IF (wndm(ji,jj) > 8.0 .AND. wndm(ji,jj) <= 25.0) THEN
430	kku=4
431	ELSE IF (wndm(ji,jj) > 25.0 ) THEN
432	kku=5
433	ENDIF
434
435	ch(ji,jj) = ( a_h(kku) + b_h(kku) * wndm(ji,jj) ** p_h(kku) &
436	+ c_h(kku) * (wndm(ji,jj)-8 ) *2) fh
437
438	ce(ji,jj) = ( a_e(kku) + b_e(kku) * wndm(ji,jj) ** p_e(kku) &
439	+ c_e(kku) * (wndm(ji,jj)-8 ) *2) fe
440
441	ch(ji,jj) = ch(ji,jj) / 1000.0
442	ce(ji,jj) = ce(ji,jj) / 1000.0
443
444	IF (wndm(ji,jj)<0.3) THEN
445	ch(ji,jj) = 1.3e-03 * fh
446	ce(ji,jj) = 1.5e-03 * fe
447	ELSE IF(wndm(ji,jj)>50.0) THEN
448	ch(ji,jj) = 1.25e-03 * fh
449	ce(ji,jj) = 1.30e-03 * fe
450	ENDIF
451
452	!!----------------------------------------------------------------------
453	!! --- calculates the SENSIBLE HEAT FLUX in MKS ( watt/m*m )
454	!!----------------------------------------------------------------------
455
456	HA(ji,jj) = rhom(ji,jj)cpch(ji,jj)wndm(ji,jj)deltemp
457
458	!!----------------------------------------------------------------------
459	!! ------------------ (iv) latent heat
460	!! --- calculates the LATENT HEAT FLUX ( watt/m*m )
461	!! --- ELAT = LrhoCe\|V\|[qs(Ts)-qa(t2d)]
462	!!----------------------------------------------------------------------
463
464	esre = shms(ji,jj) - shnow(ji,jj) ! --- calculates the term : qs(Ta)-qa(t2d)
465
466	cseep = ce(ji,jj) * wndm(ji,jj) * esre ! --- calculates the term : Ce\|V\|[qs(Ts)-qa(t2d)]
467
468	evap(ji,jj) = (cseep * rhom(ji,jj)) ! in [kg/m2/sec] !! --- calculates the EVAPORATION RATE [m/yr]
469
470	elat(ji,jj) = rhom(ji,jj) * cseep * heatlat(sst(ji,jj))
471
472	!!----------------------------------------------------------------------
473	!! --- calculates the Drag Coefficient
474	!!----------------------------------------------------------------------
475
476	!!----------------------------------------------------------------------
477	!! --- deltemp should be (Ts - Ta) in the formula estimating
478	!! --- drag coefficient
479	!!----------------------------------------------------------------------
480
481	IF( .NOT. ln_cdgw ) THEN
482	cdx(ji,jj) = cd_HR(wndm(ji,jj),deltemp)
483	ENDIF
484
485	END DO
486	END DO
487
488	!!----------------------------------------------------------------------
489	!! --- calculates the wind stresses in MKS ( newton/m*m )
490	!! --- taux= rhoCd\|V\|u tauy= rhoCd\|V\|v
491	!!----------------------------------------------------------------------
492
493	taux(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windu(:,:)
494	tauy(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windv(:,:)
495
496	CALL wrk_dealloc( jpi,jpj, rspeed, sh10now, t10now, cdx, ce, shms )
497	CALL wrk_dealloc( jpi,jpj, rhom, sstk, ch, rel_windu, rel_windv )
498	!
499	IF( nn_timing == 1 ) CALL timing_stop('fluxes_mfs')
500	!
501	END SUBROUTINE fluxes_mfs
502
503
504	REAL(wp) FUNCTION cd_HR(speed,delt)
505	!!----------------------------------------------------------------------
506	!! --- calculates the Drag Coefficient as a function of the abs. value
507	!! --- of the wind velocity ( Hellermann and Rosenstein )
508	!!----------------------------------------------------------------------
509
510	REAL(wp), INTENT(in) :: speed,delt
511	REAL(wp), PARAMETER :: a1=0.934e-3 , a2=0.788e-4, a3=0.868e-4
512	REAL(wp), PARAMETER :: a4=-0.616e-6, a5=-.120e-5, a6=-.214e-5
513
514	cd_HR = a1 + a2speed + a3delt + a4speedspeed &
515	+ a5deltdelt + a6speeddelt
516
517	END FUNCTION cd_HR
518
519	REAL(wp) function HEATLAT(t)
520	!!----------------------------------------------------------------------
521	!! --- calculates the Latent Heat of Vaporization ( J/kg ) as function of
522	!! --- the temperature ( Celsius degrees )
523	!! --- ( from A. Gill pag. 607 )
524	!!
525	!! --- Constant Latent Heat of Vaporization ( Rosati,Miyakoda 1988 )
526	!! L = 2.501e+6 (MKS)
527	!!----------------------------------------------------------------------
528
529	REAL(wp) , intent(in) :: t
530
531	heatlat = 2.5008e+6 -2.3e+3*t
532
533	END FUNCTION HEATLAT
534
535
536	SUBROUTINE qshort(hour,alat,alon,cldnow,qsw)
537	!!----------------------------------------------------------------------
538	!! * ROUTINE qshort *
539	!!
540	!! ** Purpose : Compute Solar Radiation
541	!!
542	!! ** Method : Compute Solar Radiation according Astronomical
543	!! formulae
544	!!
545	!! References : Reed RK (1975) and Reed RK (1977)
546	!!
547	!! Note: alat,alon - (lat, lon) in radians
548	!!----------------------------------------------------------------------
549	REAL(wp), INTENT (in) :: hour
550
551	REAL(wp), INTENT(in ), DIMENSION(:,:) :: alat,alon
552	REAL(wp), INTENT(in ), DIMENSION(:,:) :: cldnow
553	REAL(wp), INTENT(out), DIMENSION(:,:) :: qsw
554	REAL(wp), DIMENSION(12) :: alpham
555
556	REAL(wp), PARAMETER :: eclips=23.439* (3.141592653589793_wp / 180._wp)
557	REAL(wp), PARAMETER :: solar = 1350.
558	REAL(wp), PARAMETER :: tau = 0.7
559	REAL(wp), PARAMETER :: aozone = 0.09
560	REAL(wp), PARAMETER :: yrdays = 360.
561	REAL(wp) :: days, th0,th02,th03, sundec, thsun, coszen, qatten
562	REAL(wp) :: qzer, qdir,qdiff,qtot,tjul,sunbet
563	REAL(wp) :: albedo
564	INTEGER :: jj, ji
565
566	!!----------------------------------------------------------------------
567	!! --- albedo monthly values from Payne (1972) as means of the values
568	!! --- at 40N and 30N for the Atlantic Ocean ( hence the same latitudinal
569	!! --- band of the Mediterranean Sea ) :
570	!!----------------------------------------------------------------------
571
572	data alpham /0.095,0.08,0.065,0.065,0.06,0.06,0.06,0.06, &
573	0.065,0.075,0.09,0.10/
574
575	!!----------------------------------------------------------------------
576	!! days is the number of days elapsed until the day=nday_year
577	!!----------------------------------------------------------------------
578	days = nday_year -1.
579	th0 = 2.rpidays/yrdays
580	th02 = 2.*th0
581	th03 = 3.*th0
582
583	!! --- sun declination :
584	!!----------------------------------------------------------------------
585	sundec = 0.006918 - 0.399912cos(th0) + 0.070257sin(th0) - &
586	0.006758cos(th02) + 0.000907sin(th02) - &
587	0.002697cos(th03) + 0.001480sin(th03)
588
589	DO jj = 1, jpj
590	DO ji = 1, jpi
591
592	!! --- sun hour angle :
593	!!----------------------------------------------------------------------
594	thsun = (hour -12.)15.rad + alon(ji,jj)
595
596	!! --- cosine of the solar zenith angle :
597	!!----------------------------------------------------------------------
598	coszen =sin(alat(ji,jj))*sin(sundec) &
599	+cos(alat(ji,jj))cos(sundec)cos(thsun)
600
601	IF(coszen .LE. 5.035D-04) THEN
602	coszen = 0.0
603	qatten = 0.0
604	ELSE
605	qatten = tau**(1./coszen)
606	END IF
607
608	qzer = coszen * solar * &
609	(1.0+1.67E-2cos(rpi2.(days-3.0)/365.0))*2
610	qdir = qzer * qatten
611	qdiff = ((1.-aozone)qzer - qdir) 0.5
612	qtot = qdir + qdiff
613	tjul = (days -81.)*rad
614
615	!! --- sin of the solar noon altitude in radians :
616	!!----------------------------------------------------------------------
617	sunbet=sin(alat(ji,jj))sin(eclipssin(tjul)) + &
618	cos(alat(ji,jj))cos(eclipssin(tjul))
619
620	!! --- solar noon altitude in degrees :
621	!!----------------------------------------------------------------------
622
623	sunbet = asin(sunbet)/rad
624
625	!!----------------------------------------------------------------------
626	!! --- calculates the albedo according to Payne (1972)
627	!!----------------------------------------------------------------------
628
629	albedo = alpham(nmonth)
630
631	!!----------------------------------------------------------------------
632	!! --- ( radiation as from Reed(1977), Simpson and Paulson(1979) )
633	!! --- calculates SHORT WAVE FLUX ( watt/m*m )
634	!! --- ( Rosati,Miyakoda 1988 ; eq. 3.8 )
635	!!----------------------------------------------------------------------
636
637	IF(cldnow(ji,jj).LT.0.3) THEN
638	qsw(ji,jj) = qtot * (1.-albedo)
639	ELSE
640	qsw(ji,jj) = qtot(1.-0.62cldnow(ji,jj) &
641	+ .0019sunbet)(1.-albedo)
642	ENDIF
643
644	END DO
645	END DO
646
647	END SUBROUTINE qshort
648
649
650	!!======================================================================
651
652	END MODULE sbcblk_mfs

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

Download in other formats: