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_ecmwf.F90 in branches/2011/dev_MERCATOR_INGV_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2011/dev_MERCATOR_INGV_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_ecmwf.F90 @ 3085

Last change on this file since 3085 was 3085, checked in by cbricaud, 12 years ago
commit changes from dev_INGV_2011
File size: 30.3 KB

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

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

Download in other formats: