Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sbcblk_mfs.F90 @ 11738

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

Download in other formats: