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.

main_uvg_hpg.F90 in utils/tools/ABL_TOOLS – NEMO

Context Navigation

source: utils/tools/ABL_TOOLS/main_uvg_hpg.F90 @ 11589

Last change on this file since 11589 was 11589, checked in by gsamson, 5 years ago

dev_r11265_ABL : see #2131

ABL_TOOLS first working version (README empty and arch files ignored for now)

File size: 27.7 KB

Line
1	PROGRAM main
2	!!======================================================================
3	!! * PROGRAM main *
4	!!
5	!! ** Purpose : compute geostrophic wind or horizontal pressure gradient
6	!! from ECMWF atmospheric model vertical levels altitude,
7	!! temperature and humidity 3D fields
8	!!
9	!!======================================================================
10	!! History : 2016-10 (F. Lemarié) Original code
11	!!
12	!!----------------------------------------------------------------------
13	USE module_io ! I/O routines
14	USE module_grid ! compute input and output grids
15	!!
16	IMPLICIT NONE
17	!!----------------------------------------------------------------------
18	!!
19	!!
20	!!
21	!!----------------------------------------------------------------------
22	!
23	INTEGER :: ji,jj,jk,kt, nhym, nhyi
24	INTEGER :: jpka ! number of vertical levels for input and target grids
25	INTEGER :: jpi , jpj ! number of grid points in x and y directions
26	INTEGER :: status
27	INTEGER :: jptime,ctrl
28	INTEGER :: ioerr,ncid
29	INTEGER, ALLOCATABLE, DIMENSION(:,: ) :: ind
30	INTEGER, PARAMETER :: stdout = 6
31	!!
32	REAL(8) :: cff,tv
33	!!
34	REAL(8), ALLOCATABLE, DIMENSION(: ) :: lev ! A coefficients to reconstruct ECMWF grid
35	REAL(8), ALLOCATABLE, DIMENSION(: ) :: A_w ! A coefficients to reconstruct ECMWF grid
36	REAL(8), ALLOCATABLE, DIMENSION(: ) :: B_w ! B coefficients to reconstruct ECMWF grid
37	REAL(8), ALLOCATABLE, DIMENSION(: ) :: B_r ! B coefficients to reconstruct ECMWF grid
38	REAL(8), ALLOCATABLE, DIMENSION(: ) :: A_r ! A coefficients to reconstruct ECMWF grid
39	REAL(8), ALLOCATABLE, DIMENSION(: ) :: ph,lon,lat
40	REAL(8), ALLOCATABLE, DIMENSION(:,:,: ) :: e3t,ghw ! thickness of vertical layers in target grid
41	REAL(8), ALLOCATABLE, DIMENSION(: ) :: tmp1d, tmp_fullw, tmp_fullm ! temporary/working 1D arrays
42	REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: humi
43	REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: temp
44	REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: uhpg
45	REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: vhpg
46	REAL(8), ALLOCATABLE, DIMENSION(:,: ) :: slp, zsurf, zsurf_smth, slp_smth, ghw_smth !, slp_mask, ghw_mask
47	REAL(8), ALLOCATABLE, DIMENSION(:,: ) :: dx,dy,ff_t,tmask, tmask2
48	REAL(8), ALLOCATABLE, DIMENSION(:,: ) :: FX,FE,wrkx,wrke
49	REAL(8), ALLOCATABLE, DIMENSION(:,: ) :: dZx,dZe
50	!!
51	CHARACTER(len=500) :: file_u,file_v,file_hpg, file_geos ! ECMWF files containing wind components
52	CHARACTER(len=500) :: file_t,file_q, file_m ! ECMWF files containing tracers and mask
53	CHARACTER(len=500) :: file_z,file_p,cn_dir ! ECMWF files containing surface geopot and pressure
54	CHARACTER(len=500) :: out_file
55	CHARACTER(len=500) :: namelistf
56	CHARACTER(len=500) :: argument
57	CHARACTER(len= 20),DIMENSION(4) :: dimnames
58	CHARACTER(len= 20),DIMENSION(4) :: varnames
59	CHARACTER(6) :: mask_var ! name of mask variable in file_m file
60	!!
61	LOGICAL :: ln_read_zsurf ! read surface geopotential or not
62	LOGICAL :: ln_read_mask ! read land-sea mask or not
63	LOGICAL :: ln_perio_latbc ! use periodic BC along the domain latitudinal edges (for global data) or use zero-gradient BC (for regional data)
64	LOGICAL :: ln_c1d ! output only a single column in output file
65	LOGICAL :: ln_hpg_frc ! compute horizontal pressure gradient
66	LOGICAL :: ln_geo_wnd ! compute goestrophic wind components
67	LOGICAL :: ln_slp_smth ! apply gibbs oscillation filetring on mean sea level pressure
68	LOGICAL :: ln_drw_smth ! apply gibbs oscillation filetring on mean sea level pressure
69	LOGICAL :: ln_slp_log ! log(sea-level pressure) or sea-level pressure
70	LOGICAL :: ln_lsm_land ! if T mask is 1 over land and 0 over ocean if F it is the other way around
71	LOGICAL :: ln_impose_z1 ! impose the altitude of the first level in target grid
72	INTEGER :: ptemp_method ! way to compute potential temperature
73	! = 0 (absolute temperature)
74	! = 1 (potential temperature with local ref pressure)
75	! = 2 (potential temperature with global ref pressure on temperature perturbation)
76	! = 3 (potential temperature with global ref pressure)
77	!!
78	REAL(8), PARAMETER :: grav = 9.80665
79	REAL(8), PARAMETER :: Rd = 287.058
80	REAL(8), PARAMETER :: zvir = 0.609133
81	REAL(8), PARAMETER :: omega = 7.292116e-05
82	REAL(8), PARAMETER :: rad = 3.141592653589793 / 180.
83	REAL(8), PARAMETER :: rt = 6371229.
84	REAL(8), PARAMETER :: lat_smth = 0. !65. !!GS: possibility to smooth only above a selected latitude
85
86
87	!!---------------------------------------------------------------------
88	!! List of variables read in the namelist file
89	NAMELIST/nml_opt/ ptemp_method, ln_slp_log, ln_slp_smth, ln_read_mask, ln_perio_latbc, &
90	& ln_hpg_frc, ln_geo_wnd, ln_c1d, ln_read_zsurf, ln_lsm_land, ln_drw_smth
91	NAMELIST/nml_fld/ cn_dir, file_u, file_v, file_t, &
92	& file_q, file_z, file_p, file_hpg, file_geos, &
93	& file_m, mask_var
94	!!
95	!! get the namelist file name
96	CALL get_command_argument( 1, argument, ctrl, status)
97	!
98	SELECT CASE(status)
99	CASE(0)
100	namelistf = trim(argument)
101	CASE(-1)
102	WRITE(stdout,*) "### Error: file name too long"
103	STOP
104	CASE DEFAULT
105	namelistf = 'namelist_abl_tools'
106	END SELECT
107	!!---------------------------------------------------------------------
108
109
110	!!---------------------------------------------------------------------
111	!! read namelist variables
112	ctrl = 0
113	OPEN(50, file=namelistf, status='old', form='formatted', access='sequential', iostat=ioerr)
114	IF (ioerr /= 0) ctrl = ctrl + 1
115	READ(50,nml_opt, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
116	READ(50,nml_fld, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
117
118	IF (ctrl > 0) then
119	WRITE(stdout,*) "### E R R O R while reading namelist file '",trim(namelistf),"'"
120	WRITE(stdout,*) " ctrl = ",ctrl
121	STOP
122	ELSE
123	WRITE(stdout,*) " Namelist file ",trim(namelistf), " OK "
124	END IF
125
126	IF(ln_hpg_frc) THEN
127	out_file = trim(cn_dir)//'/'//trim(file_hpg)
128	ELSE IF(ln_geo_wnd) THEN
129	out_file = trim(cn_dir)//'/'//trim(file_geos)
130	ELSE
131	WRITE(stdout,*) "### E R R O R in namelist variable "
132	WRITE(stdout,*) "either ln_hpg_frc or ln_geo_wnd should be set to True"
133	STOP
134	END IF
135	!!---------------------------------------------------------------------
136
137
138	!!---------------------------------------------------------------------
139	! check files content
140	CALL Read_Ncdf_dim('lev',trim(cn_dir)//'/'//trim(file_t),jpka)
141	! geop_surf temp humi pressure
142	varnames = [ character(len=3) :: 'Z', 'T', 'Q', 'MSL' ]
143	IF (ln_read_zsurf) varnames(4) = "SP"
144	IF (ln_slp_log) varnames(4) = "LNSP"
145
146	ctrl = 0
147	IF (ln_read_zsurf) THEN
148	IF( .not. VAR_EXISTENCE( trim(varnames(1)) , trim(cn_dir)//'/'//trim(file_z) ) ) ctrl = ctrl + 1
149	END IF
150	IF ( .not. VAR_EXISTENCE( trim(varnames(2)) , trim(cn_dir)//'/'//trim(file_t) ) ) ctrl = ctrl + 1
151	IF ( .not. VAR_EXISTENCE( trim(varnames(3)) , trim(cn_dir)//'/'//trim(file_q) ) ) ctrl = ctrl + 1
152	IF ( .not. VAR_EXISTENCE( trim(varnames(4)) , trim(cn_dir)//'/'//trim(file_p) ) ) ctrl = ctrl + 1
153	WRITE(,) " pressure variable name: ", varnames(4)
154
155	IF ( ctrl > 0 ) THEN
156	WRITE(stdout,*) "### E R R O R while reading ECMWF atmospheric files "
157	STOP
158	ELSE
159	WRITE(stdout,*) " ECMWF atmospheric files OK "
160	END IF
161	!!---------------------------------------------------------------------
162
163
164	!!---------------------------------------------------------------------
165	!! read the dimensions for the input files
166	CALL Read_Ncdf_dim ( 'time', trim(cn_dir)//'/'//trim(file_t), jptime )
167	CALL Read_Ncdf_dim ( 'lon' , trim(cn_dir)//'/'//trim(file_t), jpi )
168	CALL Read_Ncdf_dim ( 'lat' , trim(cn_dir)//'/'//trim(file_t), jpj )
169	CALL Read_Ncdf_dim ( 'nhym' , trim(cn_dir)//'/'//trim(file_t), nhym )
170	CALL Read_Ncdf_dim ( 'nhyi' , trim(cn_dir)//'/'//trim(file_t), nhyi )
171	!
172	!!---------------------------------------------------------------------
173
174
175	!!---------------------------------------------------------------------
176	!! allocate arrays
177	ALLOCATE( A_w ( 0:jpka ) )
178	ALLOCATE( B_w ( 0:jpka ) )
179	ALLOCATE( B_r ( 1:jpka ) )
180	ALLOCATE( A_r ( 1:jpka ) )
181
182	ALLOCATE( e3t ( 1:jpi, 1:jpj, 1:jpka ) )
183	ALLOCATE( ghw ( 1:jpi, 1:jpj, 0:jpka ) )
184	ALLOCATE( slp_smth(1:jpi, 1:jpj ) )
185	ALLOCATE( ghw_smth(1:jpi, 1:jpj ) )
186	ALLOCATE( slp ( 1:jpi, 1:jpj ) )
187	ALLOCATE( zsurf ( 1:jpi, 1:jpj ) )
188	ALLOCATE( zsurf_smth ( 1:jpi, 1:jpj ) )
189	ALLOCATE( temp ( 1:jpi, 1:jpj, 1:jpka, 1) )
190	ALLOCATE( humi ( 1:jpi, 1:jpj, 1:jpka, 1) )
191	ALLOCATE( uhpg ( 1:jpi, 1:jpj, 1:jpka, 1) )
192	ALLOCATE( vhpg ( 1:jpi, 1:jpj, 1:jpka, 1) )
193	ALLOCATE( ph ( 0:jpka ) )
194	ALLOCATE( dx ( 1:jpi, 1:jpj ) )
195	ALLOCATE( dy ( 1:jpi, 1:jpj ) )
196	ALLOCATE( FX ( 1:jpi, 1:jpj ) )
197	ALLOCATE( FE ( 1:jpi, 1:jpj ) )
198	ALLOCATE( dzx ( 0:jpi, 1:jpj ) )
199	ALLOCATE( dze ( 1:jpi, 0:jpj ) )
200	ALLOCATE( wrkx ( 1:jpi, 1:jpj ) )
201	ALLOCATE( wrke ( 1:jpi, 1:jpj ) )
202	ALLOCATE( tmask ( 1:jpi, 1:jpj ) )
203	ALLOCATE( tmask2 ( 1:jpi, 1:jpj ) )
204	IF (jpka.NE.nhym) THEN
205	ALLOCATE( tmp_fullw(1:nhyi) )
206	ALLOCATE( tmp_fullm(1:nhym) )
207	END IF
208
209	!!---------------------------------------------------------------------
210	!! Read the mask and remove some closed seas
211	IF (ln_read_mask) THEN
212	CALL init_atm_mask( jpi, jpj, trim(cn_dir)//'/'//trim(file_m), trim(mask_var), ln_lsm_land, tmask)
213	ELSE
214	tmask(:,:) = 1.
215	END IF
216	tmask2 = tmask
217	!!
218
219	!! Read the static A and B coefficients for the ECMWF vertical grid
220	IF (jpka.EQ.nhym) THEN
221	CALL Read_Ncdf_var ( 'hyai', trim(cn_dir)//'/'//trim(file_t), A_w )
222	CALL Read_Ncdf_var ( 'hybi', trim(cn_dir)//'/'//trim(file_t), B_w )
223	CALL Read_Ncdf_var ( 'hyam', trim(cn_dir)//'/'//trim(file_t), A_r )
224	CALL Read_Ncdf_var ( 'hybm', trim(cn_dir)//'/'//trim(file_t), B_r )
225	ELSE
226	CALL Read_Ncdf_var ( 'hyai', trim(cn_dir)//'/'//trim(file_t), tmp_fullw )
227	A_w(0:jpka) = tmp_fullw(nhyi-(jpka+1)+1:nhyi)
228	CALL Read_Ncdf_var ( 'hybi', trim(cn_dir)//'/'//trim(file_t), tmp_fullw )
229	B_w(0:jpka) = tmp_fullw(nhyi-(jpka+1)+1:nhyi)
230	CALL Read_Ncdf_var ( 'hyam', trim(cn_dir)//'/'//trim(file_t), tmp_fullm )
231	A_r(1:jpka) = tmp_fullm(nhym-jpka+1:nhym)
232	CALL Read_Ncdf_var ( 'hybm', trim(cn_dir)//'/'//trim(file_t), tmp_fullm )
233	B_r(1:jpka) = tmp_fullm(nhym-jpka+1:nhym)
234	END IF
235
236	ALLOCATE(lat(1:jpj),lon(1:jpi),ff_t(1:jpi,1:jpj),lev(1:jpka))
237	CALL Read_Ncdf_var ( 'lon' , trim(cn_dir)//'/'//trim(file_t), lon ) !<-- longitude
238	CALL Read_Ncdf_var ( 'lat' , trim(cn_dir)//'/'//trim(file_t), lat ) !<-- latitude
239	CALL Read_Ncdf_var ( 'lev' , trim(cn_dir)//'/'//trim(file_t), lev )
240
241
242	!!---------------------------------------------------------------------
243	!++ Compute Coriolis frequency at cell centers
244	!
245	DO jj = 1, jpj
246	DO ji = 1, jpi
247	ff_t(ji,jj) = 2. * omega * SIN( rad * lat(jj) )
248	END DO
249	END DO
250
251
252	!!---------------------------------------------------------------------
253	!++ Compute dx and dy at cell centers
254	!
255	dx(:,:) = 0.
256	dy(:,:) = 0.
257	DO jj = 2, jpj-1
258	DO ji = 1, jpi-1
259	dx(ji,jj) = rt * rad * abs( lon(ji+1) - lon(ji) ) * COS( rad * lat(jj) )
260	END DO
261	END DO
262	dx( jpi,1:jpj) = dx( jpi-1,1:jpj )
263	dx(1:jpi, jpj) = dx(1:jpi , jpj-1)
264	dx(1:jpi,1 ) = dx(1:jpi ,2 )
265	!++
266	DO jj = 1, jpj-1
267	DO ji = 1, jpi
268	dy(ji,jj) = rt * rad * abs( lat(jj+1) - lat(jj) )
269	END DO
270	END DO
271	dy(1:jpi,jpj) = dy(1:jpi,jpj-1)
272
273
274	!!---------------------------------------------------------------------
275	!! create output file
276	!!
277	status = nf90_create( trim(out_file), NF90_WRITE, ncid )
278	status = nf90_close ( ncid )
279
280	CALL Write_Ncdf_dim ( 'lon' , trim(out_file), jpi )
281	CALL Write_Ncdf_dim ( 'lat' , trim(out_file), jpj )
282	CALL Write_Ncdf_dim ( 'lev' , trim(out_file), jpka )
283	CALL Write_Ncdf_dim ( 'nhym' , trim(out_file), jpka )
284	CALL Write_Ncdf_dim ( 'nhyi' , trim(out_file), jpka+1 )
285	CALL Write_Ncdf_dim ( 'time' , trim(out_file), 0 )
286
287
288	!!---------------------------------------------------------------------
289	!! Initialize the name of the dimensions for geostrophic winds in the output file
290	!!
291	dimnames(1) = 'lon'
292	dimnames(2) = 'lat'
293	dimnames(3) = 'lev'
294	dimnames(4) = 'time'
295
296	CALL Write_Ncdf_var( 'lon', 'lon', trim(out_file), lon, 'double' )
297	CALL Write_Ncdf_var( 'lat', 'lat', trim(out_file), lat, 'double' )
298	CALL Write_Ncdf_var( 'lev', 'lev', trim(out_file), lev, 'double' )
299	CALL Write_Ncdf_var( 'hyai', 'nhyi', trim(out_file), A_w, 'double' )
300	CALL Write_Ncdf_var( 'hybi', 'nhyi', trim(out_file), B_w, 'double' )
301	CALL Write_Ncdf_var( 'hyam', 'nhym', trim(out_file), A_r, 'double' )
302	CALL Write_Ncdf_var( 'hybm', 'nhym', trim(out_file), B_r, 'double' )
303
304	CALL Write_Ncdf_var( 'tmask', dimnames(1:2), trim(out_file), tmask, 'float' )
305
306
307	!!---------------------------------------------------------------------
308	! Read time variable
309	ALLOCATE(tmp1d (1:jptime))
310	CALL Read_Ncdf_var ( 'time', trim(cn_dir)//'/'//trim(file_t), tmp1d )
311	!!---------------------------------------------------------------------
312
313
314	DO kt=1,jptime
315	!
316	WRITE(stdout,*) '======================'
317	WRITE(stdout,*) 'time = ',kt,'/',jptime
318	!
319	CALL Write_Ncdf_var( 'time', dimnames(4:4), trim(out_file), tmp1d(kt:kt), kt, 'double' )
320	!
321	IF( kt == 1 ) THEN
322	CALL Duplicate_lon_lat_time( trim(cn_dir)//'/'//trim(file_t), out_file )
323	CALL Duplicate_lev_hyb ( trim(cn_dir)//'/'//trim(file_t), out_file )
324	ENDIF
325	!
326	IF ( varnames(4) == "LNSP" ) THEN
327	CALL Read_Ncdf_var( varnames(4) , trim(cn_dir)//'/'//trim(file_p), slp(:,:), kt, 1 ) !<-- log of surface pressure
328	ELSE
329	CALL Read_Ncdf_var( varnames(4) , trim(cn_dir)//'/'//trim(file_p), slp(:,:), kt )
330	END IF
331	!
332	IF (ln_slp_log) THEN
333	DO jj = 1, jpj
334	DO ji = 1, jpi
335	slp(ji,jj) = exp( slp(ji,jj) )
336	END DO
337	END DO
338	ENDIF
339	!
340	IF (ln_read_zsurf) THEN
341	CALL Read_Ncdf_var( varnames(1) , trim(cn_dir)//'/'//trim(file_z), zsurf(:,:), kt, 1 ) !<-- surface geopotential
342	ELSE
343	zsurf(:,:) = 0.
344	END IF
345	!
346	CALL Read_Ncdf_var ( varnames(2), trim(cn_dir)//'/'//trim(file_t), temp(:,:,:,:), kt ) !<-- temperature
347	CALL Read_Ncdf_var ( varnames(3), trim(cn_dir)//'/'//trim(file_q), humi(:,:,:,:), kt ) !<-- humidity
348	WHERE(humi.LT.1.E-08) humi = 1.E-08 !<-- negative values in ECMWF
349	!
350	! Smoothing of surface fields to remove gibbs oscillations (must be done on both fields or none of them)
351	!IF( ln_slp_smth ) CALL smooth_field( jpi, jpj, slp (:,:), tmask(:,:), 3 )
352	!IF (ln_read_zsurf.AND.ln_slp_smth) CALL smooth_field( jpi, jpj, zsurf(:,:), tmask(:,:), 3 )
353	!IF( ln_slp_smth ) CALL DTV_Filter( jpi, jpj, slp(:,:), tmask(:,:,1), 25, kt ) !<-- not yet robust enough
354
355	!!GS: DO NOT USE SMOOTH + LAND MASK FOR NOW (BUG)
356	IF( ln_slp_smth ) THEN
357	slp_smth(:,:) = slp(:,:)
358	wrke(:,:) = 1.
359	CALL smooth_field( jpi, jpj, slp_smth(:,1:jpj), wrke(:,1:jpj), 3 )
360	IF (ln_read_zsurf) THEN
361	zsurf_smth(:,:) = zsurf(:,:)
362	CALL smooth_field( jpi, jpj, zsurf_smth(:,:), wrke(:,:), 3 )
363	END IF
364	IF( ABS(lat_smth).GT.0.1) THEN
365	DO jj = 1, jpj
366	DO ji = 1, jpi
367	IF ((lat(jj).GE.lat_smth).OR.(lat(jj).LE.-1.*lat_smth)) THEN
368	slp(ji,jj) = slp_smth(ji,jj)
369	IF (ln_read_zsurf) zsurf(ji,jj) = zsurf_smth(ji,jj)
370	END IF
371	END DO
372	END DO
373	ELSE
374	slp(:,:) = slp_smth(:,:)
375	IF (ln_read_zsurf) zsurf(:,:) = zsurf_smth(:,:)
376	END IF
377	END IF
378	CALL Write_Ncdf_var( 'slp', dimnames(1:2), trim(out_file), slp(:,:), kt, 'float' )
379	IF (ln_read_zsurf) CALL Write_Ncdf_var( 'zsurf', dimnames(1:2), trim(out_file), zsurf(:,:), kt, 'float' )
380
381	!
382	! Compute the altitude at layer interfaces
383	ghw(:,:,1:jpka) = 0.
384	ghw(:,:, jpka) = zsurf(:,:) * (1. / grav) ! * tmask(:,:)
385	DO jj = 1, jpj
386	DO ji = 1, jpi
387	DO jk = 0, jpka
388	ph(jk) = A_w( jk ) + B_w( jk ) * slp( ji, jj ) !<-- Pa
389	END DO
390	!ph(0) = 0.1
391	IF ( nhym .EQ. jpka) ph(0) = 1.
392	DO jk = jpka,1,-1
393	tv = temp( ji, jj, jk, 1 ) * ( 1. + zvir*humi( ji, jj, jk, 1 ) ) !<-- Virtual temperature
394	e3t ( ji, jj, jk ) = (1./grav)( Rd tv * log( ph( jk ) / ph( jk-1 ) ) ) !* tmask(ji, jj)
395	ghw ( ji, jj, jk-1 ) = e3t( ji, jj, jk ) + ghw( ji, jj, jk )
396	END DO
397	END DO
398	END DO
399
400	IF( ln_slp_smth ) THEN
401	wrke(:,:) = 1.
402	DO jk = 0, jpka
403	ghw_smth(:,:) = ghw(:,:,jk)
404	CALL smooth_field( jpi, jpj, ghw_smth(:,1:jpj), wrke(:,1:jpj), 3 )
405	IF( ABS(lat_smth).GT.0.1) THEN
406	DO jj = 1, jpj
407	DO ji = 1, jpi
408	IF ((lat(jj).GE.lat_smth).OR.(lat(jj).LE.-1.*lat_smth)) ghw(ji,jj,jk) = ghw_smth(ji,jj)
409	END DO
410	END DO
411	ELSE
412	ghw(:,:,jk) = ghw_smth(:,:)
413	END IF
414	END DO
415	END IF
416
417	!
418	! Compute horizontal gradient of slp in x-direction (FX = dslp / dx)
419	FX(:,:) = 0.
420	DO jj = 1, jpj
421	DO ji = 1, jpi-1
422	IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji+1,jj) .gt. 0.5)) THEN
423	cff = 2. / ( dx( ji+1,jj ) + dx( ji,jj ) )
424	FX( ji ,jj ) = cff * ( slp( ji+1,jj ) - slp( ji,jj ) )
425	ELSE
426	tmask2(ji:ji+1,jj) = 0.
427	END IF
428	END DO
429	END DO
430
431	IF (ln_perio_latbc) THEN
432	! apply periodicity
433	DO jj = 1, jpj
434	IF ((tmask(1,jj) .gt. 0.5).AND.(tmask(jpi,jj) .gt. 0.5)) THEN
435	cff = 2. / ( dx( 1, jj) + dx( jpi, jj) )
436	FX( jpi , jj) = cff * ( slp( 1, jj) - slp( jpi, jj) )
437	ELSE
438	tmask2( 1, jj) = 0.
439	tmask2( jpi, jj) = 0.
440	END IF
441	END DO
442	ELSE
443	! apply no-gradient
444	DO jj = 1, jpj
445	FX( jpi ,jj ) = FX( jpi-1 ,jj )
446	END DO
447	ENDIF
448
449	!
450	! Compute horizontal gradient of slp in y-direction (FE = dslp / dy)
451	FE(:,:) = 0.
452	DO jj = 1, jpj-1
453	DO ji = 1, jpi
454	IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji,jj+1) .gt. 0.5)) THEN
455	cff = 2. / ( dy( ji,jj+1 ) + dy( ji,jj ) )
456	FE( ji ,jj ) = cff * ( slp( ji,jj+1 ) - slp( ji,jj ) )
457	ELSE
458	tmask2(ji,jj:jj+1) = 0.
459	END IF
460	END DO
461	END DO
462
463	! apply no-gradient
464	DO ji = 1, jpi
465	FE( ji ,jpj ) = FE( ji ,jpj-1 )
466	END DO
467
468	!
469	!++ Compute the geostrophic winds
470	!
471	dZx(:,:) = 0.
472	dZe(:,:) = 0.
473	wrkX(:,:) = 0.
474	wrkE(:,:) = 0.
475	!////////////
476	DO jk=1,jpka
477	!////////////
478
479	!
480	! Compute horizontal gradient of altitude in x-direction along the coordinate dZx = (dz / dx)s
481	DO jj = 1, jpj
482	DO ji = 1, jpi-1
483	IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji+1,jj) .gt. 0.5)) THEN
484	cff = 1. / (dx(ji+1,jj)+dx(ji,jj))
485	dZx(ji,jj) = cff * ( (ghw( ji+1, jj, jk-1) - ghw( ji, jj, jk-1)) &
486	& +(ghw( ji+1, jj, jk ) - ghw( ji, jj, jk )) )
487	ELSE
488	tmask2(ji:ji+1,jj) = 0.
489	END IF
490	END DO
491	END DO
492
493	IF (ln_perio_latbc) THEN
494	! apply periodicity
495	DO jj = 1, jpj
496	IF ((tmask(1,jj) .gt. 0.5).AND.(tmask(jpi,jj) .gt. 0.5)) THEN
497	cff = 1. / (dx(1,jj)+dx(jpi,jj))
498	dZx(jpi,jj) = cff * ( (ghw( 1, jj, jk-1) - ghw( jpi, jj, jk-1)) &
499	& +(ghw( 1, jj, jk ) - ghw( jpi, jj, jk )) )
500	dZx( 0,jj) = dZx(jpi,jj)
501	ELSE
502	tmask2( 1, jj) = 0.
503	tmask2(jpi, jj) = 0.
504	END IF
505	END DO
506	ELSE
507	! apply no-gradient
508	DO jj = 1, jpj
509	dZx(jpi,jj) = dZx(jpi-1,jj)
510	dZx( 0,jj) = dZx( 1,jj)
511	END DO
512	END IF
513
514	!
515	! Compute horizontal gradient of altitude in y-direction along the coordinate dZy = (dz / dy)s
516	DO jj = 1, jpj-1
517	DO ji = 1, jpi
518	IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji,jj+1) .gt. 0.5)) THEN
519	cff = 1. / (dy(ji,jj)+dy(ji,jj+1))
520	dZe(ji,jj) = cff * ( (ghw( ji, jj+1, jk-1) - ghw( ji, jj, jk-1)) &
521	& +(ghw( ji, jj+1, jk ) - ghw( ji, jj, jk )) )
522	ELSE
523	tmask2(ji,jj:jj+1) = 0.
524	END IF
525	END DO
526	END DO
527
528	! apply no-gradient
529	DO ji = 1, jpi
530	dZe(ji,jpj) = dZe(ji,jpj-1)
531	dZe(ji, 0) = dZe(ji, 1)
532	END DO
533
534
535	! Compute horizontal pressure gradient in x-direction along the coordinate wrkX = (dp/dx)s
536	DO jj = 1, jpj
537	DO ji = 2, jpi
538	IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji-1,jj) .gt. 0.5)) THEN
539	cff = slp(ji,jj) * (B_w(jk)-B_w(jk-1)) + (A_w(jk)-A_w(jk-1))
540	wrkX(ji,jj) = B_r(jk) * 0.5 * (FX(ji,jj)+FX(ji-1,jj)) &
541	& * (ghw(ji,jj,jk)-ghw(ji,jj,jk-1)) / cff
542	ELSE
543	tmask2(ji-1:ji,jj) = 0.
544	END IF
545	END DO
546	END DO
547
548	IF (ln_perio_latbc) THEN
549	! apply periodicity
550	ji = 1
551	DO jj = 1, jpj
552	IF ((tmask(1,jj) .gt. 0.5).AND.(tmask(jpi,jj) .gt. 0.5)) THEN
553	cff = slp(ji,jj) * (B_w(jk)-B_w(jk-1)) + (A_w(jk)-A_w(jk-1))
554	wrkX(ji,jj) = B_r(jk) * 0.5 * (FX(ji,jj)+FX(jpi,jj)) &
555	& * (ghw(ji,jj,jk)-ghw(ji,jj,jk-1)) / cff
556	ELSE
557	tmask2( 1,jj) = 0.
558	tmask2(jpi,jj) = 0.
559	END IF
560	END DO
561	ELSE
562	! apply no gradient
563	DO jj = 1, jpj
564	wrkX(1,jj) = wrkX(2,jj)
565	END DO
566	END IF
567
568
569	! Compute horizontal pressure gradient in y-direction along the coordinate wrkE = (dp/dy)s
570	DO jj = 2, jpj
571	DO ji = 1, jpi
572	IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji,jj-1) .gt. 0.5)) THEN
573	cff = slp(ji,jj) * (B_w(jk)-B_w(jk-1)) + (A_w(jk)-A_w(jk-1))
574	wrkE(ji,jj) = B_r(jk) * 0.5 * (FE(ji,jj)+FE(ji,jj-1)) &
575	& * (ghw(ji,jj,jk)-ghw(ji,jj,jk-1)) / cff
576	ELSE
577	tmask2(ji,jj-1:jj) = 0.
578	END IF
579	END DO
580	END DO
581
582	! apply no gradient
583	jj = 1
584	DO ji = 1, jpi
585	wrkE(ji,1) = wrkE(ji,2)
586	END DO
587
588
589	!+++ Finalize pressure gradient/geostrophic wind computation
590	IF(ln_hpg_frc) THEN
591	DO jj=1,jpj
592	DO ji=1,jpi
593	IF (tmask2(ji,jj).GT.0.5) THEN
594	uhpg(ji,jj,jk,1) = - grav( wrkX(ji,jj) - 0.5( dZx(ji,jj)+dZx(ji-1,jj ) ) )
595	if (lat(1).GT.0.) vhpg(ji,jj,jk,1) = grav( wrkE(ji,jj) - 0.5( dZe(ji,jj)+dZe(ji ,jj-1) ) )
596	if (lat(1).LT.0.) vhpg(ji,jj,jk,1) = - grav( wrkE(ji,jj) - 0.5( dZe(ji,jj)+dZe(ji ,jj-1) ) )
597	ELSE
598	uhpg(ji,jj,jk,1) = 0.
599	vhpg(ji,jj,jk,1) = 0.
600	END IF
601	END DO
602	END DO
603	ELSE
604	DO jj=1,jpj
605	DO ji=1,jpi
606	IF (tmask2(ji,jj).GT.0.5) THEN
607	cff = grav / ff_t(ji,jj)
608	! geostrophic wind computed only where Coriolis .ge. 3.e-5 (~12deg)
609	if(abs(ff_t(ji,jj)) < 2.5e-5) cff = 0.
610	! minus sign for uhpg because y-derivatives are inverted
611	vhpg(ji,jj,jk,1) = - cff( wrkX(ji,jj) - 0.5( dZx(ji,jj)+dZx(ji-1,jj ) ) )
612	if (lat(1).GT.0.) uhpg(ji,jj,jk,1) = - cff( wrkE(ji,jj) - 0.5( dZe(ji,jj)+dZe(ji ,jj-1) ) )
613	if (lat(1).LT.0.) uhpg(ji,jj,jk,1) = cff( wrkE(ji,jj) - 0.5( dZe(ji,jj)+dZe(ji ,jj-1) ) )
614	ELSE
615	uhpg(ji,jj,jk,1) = 0.
616	vhpg(ji,jj,jk,1) = 0.
617	END IF
618	END DO
619	END DO
620	ENDIF
621
622	!////////////
623	END DO ! jk
624	!////////////
625
626	IF (ln_geo_wnd) THEN
627	CALL Write_Ncdf_var ( 'ugeo', dimnames(1:4), trim(out_file), uhpg, kt, 'float' )
628	CALL Write_Ncdf_var ( 'vgeo', dimnames(1:4), trim(out_file), vhpg, kt, 'float' )
629	ELSE
630	CALL Write_Ncdf_var ( 'uhpg', dimnames(1:4), trim(out_file), uhpg, kt, 'float' )
631	CALL Write_Ncdf_var ( 'vhpg', dimnames(1:4), trim(out_file), vhpg, kt, 'float' )
632	END IF
633
634	IF (kt .EQ. 1) THEN
635	tmask2(:,1) = 0. ! force northern v line drowning
636	tmask2(:,jpj) = 0. ! force northern v line drowning
637	CALL Write_Ncdf_var( 'tmask', dimnames(1:2), trim(out_file), tmask2, 'float' )
638	END IF
639
640	END DO ! kt
641	!
642	DEALLOCATE( zsurf, zsurf_smth, slp, slp_smth, temp,humi,uhpg,vhpg)
643	IF (jpka.NE.nhym) DEALLOCATE(tmp_fullw,tmp_fullm)
644	!
645	STOP
646	!
647	END PROGRAM main

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

Download in other formats: