Context Navigation

sbcblk_mfs.F90 @ 9383

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

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

Download in other formats: