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/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90 @ 9366

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

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

Download in other formats: