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_vinterp.F90 in utils/tools_ticket2457/ABL_TOOLS – NEMO

Context Navigation

source: utils/tools_ticket2457/ABL_TOOLS/main_vinterp.F90 @ 14431

Last change on this file since 14431 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.6 KB

Line
1	PROGRAM main
2	!!======================================================================
3	!! * PROGRAM main *
4	!!
5	!! ** Purpose : Vertical interpolation of ECMWF dataset on a given fixed
6	!! vertical grid
7	!!
8	!!======================================================================
9	!! History : 2016-10 (F. Lemarié) Original code
10	!!
11	!!----------------------------------------------------------------------
12	USE module_io ! I/O routines
13	USE module_interp ! vertical interpolation 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, jk_in, nhym, nhyi
24	INTEGER :: jpka_in, 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 :: iloc, jloc ! grid indexes for c1d case
27	INTEGER :: status
28	INTEGER :: jptime,ctrl
29	INTEGER :: ioerr
30	INTEGER, ALLOCATABLE, DIMENSION(:,: ) :: ind
31	INTEGER, PARAMETER :: stdout = 6
32	INTEGER, PARAMETER :: jp_weno = 1
33	INTEGER, PARAMETER :: jp_spln = 2
34	!!
35	REAL(8) :: hc,hmax,theta_s,z1 ! parameters related to the target vertical grid
36	REAL(8) :: cff
37	!!
38	REAL(8), ALLOCATABLE, DIMENSION(: ) :: A_w ! A coefficients to reconstruct ECMWF grid
39	REAL(8), ALLOCATABLE, DIMENSION(: ) :: A_wa ! A coefficients to reconstruct ECMWF grid
40	REAL(8), ALLOCATABLE, DIMENSION(: ) :: B_w ! B coefficients to reconstruct ECMWF grid
41	REAL(8), ALLOCATABLE, DIMENSION(: ) :: B_wa ! B coefficients to reconstruct ECMWF grid
42	REAL(8), ALLOCATABLE, DIMENSION(: ) :: tmp1d, tmp_fullw, tmp_fullm ! temporary/working 1D arrays
43	REAL(8), ALLOCATABLE, DIMENSION(: ) :: e3t,e3w ! thickness of vertical layers in target grid
44	REAL(8), ALLOCATABLE, DIMENSION(: ) :: ght,ghw ! altitude of vertical grid points
45	REAL(8), ALLOCATABLE, DIMENSION(:,: ) :: e3_bak
46	REAL(8), ALLOCATABLE, DIMENSION(:,:,: ) :: ghw_in ! altitude of cell interfaces of ECMWF grid
47	REAL(8), ALLOCATABLE, DIMENSION(:,:,: ) :: e3t_in ! thickness of vertical layers in ECMWF grid
48	REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: humi
49	REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: tair, tpot
50	REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: varout, varc1d
51	REAL(8), ALLOCATABLE, DIMENSION(:,: ) :: slp, zsurf
52	REAL(8), ALLOCATABLE, DIMENSION(:,: ) :: tmask, tmask2 ! land-sea mask
53	!!
54	CHARACTER(len=500) :: file_u,file_v,file_hpg,file_geos ! ECMWF files containing wind components
55	CHARACTER(len=500) :: file_t,file_q, file_m ! ECMWF files containing tracers and mask
56	CHARACTER(len=500) :: file_z,file_p,cn_dir,file_in ! ECMWF files containing surface geopot and pressure
57	CHARACTER(len=500) :: grd_file, abl_file, drwn_file, out_file
58	CHARACTER(len=500) :: namelistf,stmp
59	CHARACTER(len=500) :: argument, var_file
60	CHARACTER(len= 20),DIMENSION(4) :: dimnames
61	CHARACTER(len= 20),DIMENSION(11) :: varnames, outnames
62	CHARACTER(len=500),DIMENSION(11) :: filnames
63	CHARACTER(6) :: mask_var ! name of mask variable in file_m file
64	CHARACTER(6) :: var_name
65	!!
66	LOGICAL :: ln_read_zsurf ! read surface geopotential or not
67	LOGICAL :: ln_read_mask ! read land-sea mask or not
68	LOGICAL :: ln_perio_latbc ! use periodic BC along the domain latitudinal edges (for global data) or use zero-gradient BC (for regional data)
69	LOGICAL :: ln_c1d ! output only a single column in output file
70	LOGICAL :: ln_hpg_frc ! compute horizontal pressure gradient
71	LOGICAL :: ln_geo_wnd ! compute goestrophic wind components
72	LOGICAL :: ln_slp_smth ! apply gibbs oscillation filetring on mean sea level pressure
73	LOGICAL :: ln_drw_smth ! apply gibbs oscillation filetring on mean sea level pressure
74	LOGICAL :: ln_slp_log ! log(sea-level pressure) or sea-level pressure
75	LOGICAL :: ln_lsm_land ! if T mask is 1 over land and 0 over ocean if F it is the other way around
76	LOGICAL :: ln_impose_z1 ! impose the altitude of the first level in target grid
77	INTEGER :: ptemp_method ! way to compute potential temperature
78	! = 0 (absolute temperature)
79	! = 1 (potential temperature with local ref pressure)
80	! = 2 (potential temperature with global ref pressure on temperature perturbation)
81	! = 3 (potential temperature with global ref pressure)
82	!!
83	REAL(8), PARAMETER :: grav = 9.80665
84
85	!!---------------------------------------------------------------------
86	!! List of variables read in the namelist file
87	NAMELIST/nml_dom/ jpka, hmax, theta_s, hc, ln_impose_z1, z1
88	NAMELIST/nml_opt/ ptemp_method, ln_slp_log, ln_slp_smth, ln_read_mask, ln_perio_latbc, &
89	& ln_hpg_frc, ln_geo_wnd, ln_c1d, ln_read_zsurf, ln_lsm_land, ln_drw_smth
90	NAMELIST/nml_fld/ cn_dir, file_u, file_v, file_t, &
91	& file_q, file_z, file_p, file_hpg, file_geos, &
92	& file_m, mask_var
93	NAMELIST/nml_out/ grd_file, abl_file, drwn_file, var_name
94	NAMELIST/nml_c1d/ iloc, jloc
95	!!
96	!! get the namelist file name
97	CALL get_command_argument( 1, argument, ctrl, status)
98	!
99	SELECT CASE(status)
100	CASE(0)
101	namelistf = trim(argument)
102	CASE(-1)
103	WRITE(stdout,*) "### Error: file name too long"
104	STOP
105	CASE DEFAULT
106	namelistf = 'namelist_abl_tools'
107	END SELECT
108	!!---------------------------------------------------------------------
109
110
111	!!---------------------------------------------------------------------
112	!! read namelist variables
113	ctrl = 0
114	OPEN(50, file=namelistf, status='old', form='formatted', access='sequential', iostat=ioerr)
115	IF (ioerr /= 0) ctrl = ctrl + 1
116	READ(50,nml_dom, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
117	READ(50,nml_opt, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
118	READ(50,nml_fld, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
119	READ(50,nml_out, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
120	IF( ln_c1d ) THEN
121	print*,'c1d is activated'
122	READ(50,nml_c1d, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
123	ENDIF
124
125	IF (ctrl > 0) THEN
126	WRITE(stdout,*) "### E R R O R while reading namelist file '",trim(namelistf),"'"
127	WRITE(stdout,*) " ctrl = ",ctrl
128	STOP
129	ELSE
130	WRITE(stdout,*) " Namelist file ",trim(namelistf)," OK "
131	ENDIF
132	IF( ln_hpg_frc .AND. ln_geo_wnd ) THEN
133	WRITE(stdout,*) "### E R R O R conflicting options "
134	WRITE(stdout,*) "ln_hpg_frc and ln_geo_wnd can not both be set to True"
135	STOP
136	ENDIF
137
138	SELECT CASE (ptemp_method)
139	CASE(0)
140	WRITE(stdout,*) "Absolute temperature option is activated"
141	CASE(1)
142	WRITE(stdout,*) "Potential temperature option with local reference pressure is activated"
143	CASE(2)
144	WRITE(stdout,*) "Potential temperature option with global reference pressure on temperature perturbation is activated"
145	CASE(3)
146	WRITE(stdout,*) "Potential temperature option with global reference pressure is activated"
147	END SELECT
148
149	IF(ln_slp_smth) WRITE(stdout,*) "MSLP smoothing option is activated"
150	IF(ln_hpg_frc ) WRITE(stdout,*) "Large-scale pressure gradient will be interpolated"
151	IF(ln_geo_wnd ) WRITE(stdout,*) "Geostrophic winds will be interpolated"
152	!!---------------------------------------------------------------------
153
154	!!-------------------------------------------------------------------------------------
155	!! list of variables to treat
156	!!
157	!! get the variable name
158	CALL get_command_argument( 2, argument, ctrl, status)
159	SELECT CASE(status)
160	CASE(0)
161	var_name = trim(argument)
162	CASE(-1)
163	WRITE(stdout,*) "### Error: file name too long"
164	STOP
165	END SELECT
166	WRITE(stdout,*) "var_name: ", trim(var_name), "lenght: ", Len_Trim(var_name)
167	!!
168	IF( ln_hpg_frc ) THEN
169	file_in = trim(file_hpg )
170	ELSE
171	file_in = trim(file_geos)
172	ENDIF
173	varnames = [character(len=4 ) :: 'Z' , 'T', 'Q', 'U', 'V', 'MSL', 'LSM', 'uhpg', 'vhpg', 'ugeo', 'vgeo' ]
174	outnames = [character(len=5 ) :: 'zsurf', 'tair', 'humi', 'uwnd', 'vwnd', 'slp', '', 'uhpg', 'vhpg', 'ugeo', 'vgeo' ]
175	filnames = [character(len=500) :: file_z, file_t, file_q, file_u, file_v, file_p, file_m, file_in, file_in, file_in, file_in ]
176	IF( ln_slp_log ) varnames(6) = 'LNSP'
177
178	!!---------------------------------------------------------------------
179	! check files content
180	ctrl = 0
181	CALL Read_Ncdf_dim('lev',trim(cn_dir)//'/'//trim(file_t),jpka_in)
182	!
183	IF (ln_read_zsurf) THEN
184	IF( .not. VAR_EXISTENCE( trim(varnames(1)) , trim(cn_dir)//'/'//trim(file_z) ) &
185	& .or. jpka_in == 1 ) ctrl = ctrl + 1
186	ENDIF
187	WRITE(stdout,*) trim(varnames(1)), ctrl
188	IF ( .not. VAR_EXISTENCE( trim(varnames(2)) , trim(cn_dir)//'/'//trim(file_t) ) ) ctrl = ctrl + 1
189	WRITE(stdout,*) trim(varnames(2)), ctrl
190	IF ( .not. VAR_EXISTENCE( trim(varnames(3)) , trim(cn_dir)//'/'//trim(file_q) ) ) ctrl = ctrl + 1
191	WRITE(stdout,*) trim(varnames(3)), ctrl
192	IF ( .not. VAR_EXISTENCE( trim(varnames(4)) , trim(cn_dir)//'/'//trim(file_u) ) ) ctrl = ctrl + 1
193	WRITE(stdout,*) trim(varnames(4)), ctrl
194	IF ( .not. VAR_EXISTENCE( trim(varnames(5)) , trim(cn_dir)//'/'//trim(file_v) ) ) ctrl = ctrl + 1
195	WRITE(stdout,*) trim(varnames(5)), ctrl
196	IF ( .not. VAR_EXISTENCE( trim(varnames(6)) , trim(cn_dir)//'/'//trim(file_p) ) ) ctrl = ctrl + 1
197	WRITE(stdout,*) trim(varnames(6)), ctrl
198	IF (ln_read_mask) THEN
199	IF ( .not. VAR_EXISTENCE( trim(varnames(7)) , trim(cn_dir)//'/'//trim(file_m) ) ) ctrl = ctrl + 1
200	WRITE(stdout,*) trim(varnames(7)), ctrl
201	varnames(7) = TRIM(mask_var)
202	ENDIF
203	IF (ln_hpg_frc) THEN
204	IF ( .not. VAR_EXISTENCE( trim(varnames(8)) , trim(cn_dir)//'/'//trim(file_hpg) ) ) ctrl = ctrl + 1
205	WRITE(stdout,*) trim(varnames(8)), ctrl
206	IF ( .not. VAR_EXISTENCE( trim(varnames(9)) , trim(cn_dir)//'/'//trim(file_hpg) ) ) ctrl = ctrl + 1
207	WRITE(stdout,*) trim(varnames(9)), ctrl
208	ENDIF
209	IF (ln_geo_wnd) THEN
210	IF ( .not. VAR_EXISTENCE( trim(varnames(10)) , trim(cn_dir)//'/'//trim(file_geos) ) ) ctrl = ctrl + 1
211	WRITE(stdout,*) trim(varnames(10)), ctrl
212	IF ( .not. VAR_EXISTENCE( trim(varnames(11)) , trim(cn_dir)//'/'//trim(file_geos) ) ) ctrl = ctrl + 1
213	WRITE(stdout,*) trim(varnames(11)), ctrl
214	ENDIF
215
216	IF ( ctrl > 0 ) THEN
217	WRITE(stdout,*) "### E R R O R while reading ECMWF atmospheric files "
218	STOP
219	ELSE
220	WRITE(stdout,*) " ECMWF atmospheric files OK "
221	ENDIF
222	!!---------------------------------------------------------------------
223
224
225	!!---------------------------------------------------------------------
226	!! read the dimensions for the input files
227	CALL Read_Ncdf_dim ( 'time', trim(cn_dir)//'/'//trim(file_t), jptime )
228	CALL Read_Ncdf_dim ( 'lon' , trim(cn_dir)//'/'//trim(file_t), jpi )
229	CALL Read_Ncdf_dim ( 'lat' , trim(cn_dir)//'/'//trim(file_t), jpj )
230	CALL Read_Ncdf_dim ( 'nhym' , trim(cn_dir)//'/'//trim(file_t), nhym )
231	CALL Read_Ncdf_dim ( 'nhyi' , trim(cn_dir)//'/'//trim(file_t), nhyi )
232	!WRITE(stdout,*) "jpka_in, jptime, jpi, jpj, nhym, nhyi: ", jpka_in, jptime, jpi, jpj, nhym, nhyi
233	!
234	!!---------------------------------------------------------------------
235
236
237	!!---------------------------------------------------------------------
238	!! allocate arrays
239	ALLOCATE( A_w ( 0:jpka_in) )
240	ALLOCATE( B_w ( 0:jpka_in) )
241	ALLOCATE( e3t_in ( 1:jpi, 1:jpj, 1:jpka_in) )
242	ALLOCATE( ghw_in ( 1:jpi, 1:jpj, 0:jpka_in) )
243	ALLOCATE( slp ( 1:jpi, 1:jpj ) )
244	ALLOCATE( zsurf ( 1:jpi, 1:jpj ) )
245	ALLOCATE( tair ( 1:jpi, 1:jpj, 1:jpka_in, 1) )
246	ALLOCATE( tpot ( 1:jpi, 1:jpj, 1:jpka_in, 1) )
247	ALLOCATE( humi ( 1:jpi, 1:jpj, 1:jpka_in, 1) )
248	ALLOCATE( varout ( 1:jpi, 1:jpj, 1:jpka+1 , 1) )
249	ALLOCATE( ind ( 1:jpi, 1:jpj ) )
250	ALLOCATE( e3_bak ( 1:jpi, 1:jpj ) )
251	ALLOCATE( ght ( 1:jpka+1 ) )
252	ALLOCATE( ghw ( 1:jpka+1 ) )
253	ALLOCATE( e3t ( 1:jpka+1 ) )
254	ALLOCATE( e3w ( 1:jpka+1 ) )
255	ALLOCATE( tmask ( 1:jpi, 1:jpj ) )
256	ALLOCATE( tmask2 ( 1:jpi, 1:jpj ) )
257	IF( ln_c1d ) ALLOCATE( varc1d( 1:3, 1:3, 1:jpka+1 , 1 ) )
258	IF( ln_c1d ) varc1d( 1:3, 1:3, 1:jpka+1 , 1 ) = 0.
259	IF (jpka_in.NE.nhym) THEN
260	ALLOCATE( tmp_fullw(1:nhyi) )
261	ALLOCATE( tmp_fullm(1:nhym) )
262	ENDIF
263	!
264	varout(:,:,:,1) = 0.
265
266	!!---------------------------------------------------------------------
267	!! Read the mask and remove some closed seas
268	IF (ln_read_mask) THEN
269	CALL init_atm_mask(jpi,jpj,trim(cn_dir)//'/'//trim(file_m),trim(mask_var),ln_lsm_land,tmask)
270	IF (ln_geo_wnd) CALL Read_Ncdf_var ( 'tmask' , trim(cn_dir)//'/'//trim(file_geos), tmask2(:,:) )
271	IF (ln_hpg_frc) CALL Read_Ncdf_var ( 'tmask' , trim(cn_dir)//'/'//trim(file_hpg) , tmask2(:,:) )
272	ELSE
273	tmask(:,:) = 1.
274	tmask2(:,:) = 1.
275	ENDIF
276	!!
277
278	!!---------------------------------------------------------------------
279	!! Compute the altitude and layer thickness of the target grid
280	CALL init_target_grid ( jpka, ght, ghw, e3t, e3w, hmax, hc, theta_s, &
281	& ln_impose_z1, z1 )
282
283	!! Write the grid file for the target grid
284	CALL Write_Grid_File ( jpka, ght, ghw, e3t, e3w, trim(cn_dir)//'/'//trim(grd_file) )
285
286
287	!! Read the static A and B coefficients for the ECMWF vertical grid
288	IF (jpka_in.EQ.nhym) THEN
289	CALL Read_Ncdf_var ( 'hyai', trim(cn_dir)//'/'//trim(file_t), A_w )
290	CALL Read_Ncdf_var ( 'hybi', trim(cn_dir)//'/'//trim(file_t), B_w )
291	ELSE
292	CALL Read_Ncdf_var ( 'hyai', trim(cn_dir)//'/'//trim(file_t), tmp_fullw )
293	A_w(0:jpka_in) = tmp_fullw(nhyi-(jpka_in+1)+1:nhyi)
294	CALL Read_Ncdf_var ( 'hybi', trim(cn_dir)//'/'//trim(file_t), tmp_fullw )
295	B_w(0:jpka_in) = tmp_fullw(nhyi-(jpka_in+1)+1:nhyi)
296	ENDIF
297
298
299	!!---------------------------------------------------------------------
300	!! create output file
301	!!
302	IF (Len_Trim(var_name) == 0) THEN
303	out_file = trim(cn_dir)//'/'//trim(abl_file)
304	ELSE
305	out_file = trim(cn_dir)//'/'//trim(var_name)//'_'//trim(abl_file)
306	ENDIF
307
308	IF(ln_c1d) THEN
309	CALL Init_output_File_c1d ( jpi, jpj, jpka, trim(cn_dir)//'/'//trim(file_t), out_file, tmask(:,:), iloc, jloc )
310	ELSE
311	CALL Init_output_File ( jpi, jpj, jpka, trim(cn_dir)//'/'//trim(file_t), out_file, tmask(:,:) )
312	ENDIF
313
314	!!---------------------------------------------------------------------
315	!! Initialize the name of the dimensions for the result of the interpolation
316	!!
317	dimnames(1) = 'lon'
318	dimnames(2) = 'lat'
319	dimnames(3) = 'jpka'
320	dimnames(4) = 'time'
321	CALL Write_Ncdf_var( 'tmask', dimnames(1:2), trim(out_file), tmask2, 'float' )
322
323	!!---------------------------------------------------------------------
324	! Read time variable
325	ALLOCATE(tmp1d (1:jptime))
326	CALL Read_Ncdf_var ( 'time', trim(cn_dir)//'/'//trim(file_t), tmp1d )
327	!!---------------------------------------------------------------------
328
329	DO kt=1,jptime
330	!
331	WRITE(stdout,*) '======================'
332	WRITE(stdout,*) 'time = ',kt,'/',jptime
333	!
334	CALL Write_Ncdf_var( 'time', dimnames(4:4), trim(out_file), tmp1d(kt:kt), kt, 'double' )
335	!
336	IF( kt == 1 ) THEN
337	CALL Duplicate_lon_lat_time( trim(cn_dir)//'/'//trim(file_t), out_file )
338	CALL add_globatt_real( out_file, "jpka" , REAL(jpka) )
339	CALL add_globatt_real( out_file, "hmax" , hmax )
340	CALL add_globatt_real( out_file, "theta_s", theta_s )
341	CALL add_globatt_real( out_file, "hc" , hc )
342	IF (ln_impose_z1) CALL add_globatt_real( out_file, "z1", z1 )
343	ENDIF
344	!
345	! read SLP
346	!CALL Read_Ncdf_var ( varnames(6) , trim(cn_dir)//'/'//trim(file_p), slp, kt ) !<-- (log of) surface pressure
347	!IF (ln_slp_log) THEN
348	! DO jj = 1, jpj
349	! DO ji = 1, jpi
350	! slp( ji ,jj ) = exp( slp(ji ,jj ) )
351	! END DO
352	! END DO
353	!ENDIF
354	CALL Read_Ncdf_var ( outnames(6) , trim(cn_dir)//'/'//trim(file_in), slp, kt ) !<-- (log of) surface pressure
355	!
356	! read ZSURF
357	IF (ln_read_zsurf) THEN
358	!CALL Read_Ncdf_var ( varnames(1) , trim(cn_dir)//'/'//trim(file_z), zsurf, kt ) !<-- surface geopotential
359	CALL Read_Ncdf_var ( outnames(1) , trim(cn_dir)//'/'//trim(file_in), zsurf, kt ) !<-- surface geopotential
360	ELSE
361	zsurf(:,:) = 0.
362	ENDIF
363	!
364	! Smoothing of SLP and ZSURF to remove gibbs oscillations (must be done on both fields or none of them)
365	!IF( ln_slp_smth ) CALL smooth_field( jpi, jpj, slp(:,:), tmask(:,:), 3 )
366	!IF (ln_read_zsurf.AND.ln_slp_smth) CALL smooth_field( jpi, jpj, zsurf(:,:), tmask(:,:), 3 )
367	!IF( ln_slp_smth ) CALL DTV_Filter( jpi, jpj, slp(:,:), tmask(:,:), 25, kt ) !<-- not yet robust enough
368	!
369	! read tair and HUMI
370	CALL Read_Ncdf_var ( varnames(2) , trim(cn_dir)//'/'//trim(file_t), tair, kt ) !<-- temperature
371	CALL Read_Ncdf_var ( varnames(3) , trim(cn_dir)//'/'//trim(file_q), humi, kt ) !<-- humidity
372	WHERE(humi.LT.1.E-08) varout = 1.E-08 !<-- negative values in ECMWF
373	!
374	! Reconstruct the ERA-Interim vertical grid in terms of altitude
375	ghw_in(:,:,1:jpka_in) = 0.
376	ghw_in(:,:, jpka_in) = zsurf(:,:) * (1. / grav)
377	CALL get_atm_grid( jpi, jpj, jpka_in, slp, tair, &
378	& humi, A_w, B_w, e3t_in, ghw_in )
379	!
380	! Compute potential temperature
381	tpot = tair ! save tpot
382	CALL get_pot_temp ( jpi, jpj, jpka_in, slp, tpot, A_w, B_w, &
383	& tmask(:,:), ptemp_method, humi(:,:,jpka_in,1), 0.5*ghw_in(:,:,jpka_in-1) )
384
385	! Flip the vertical axis to go from k=0 at the bottom to k=N_in at the top of the atmosphere
386	CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, e3t_in )
387	CALL flip_vert_dim ( 0, jpka_in, jpi, jpj, ghw_in )
388	CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, tpot(:,:,:,1) )
389	CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, humi(:,:,:,1) )
390
391	! Correct the layer thickness to match hmax
392	DO jj = 1, jpj
393	DO ji = 1, jpi
394	cff = 0.
395	DO jk=1,jpka_in
396	cff = cff + e3t_in( ji, jj, jk )
397	IF ( cff > hmax ) THEN
398	jk_in = jk
399	EXIT
400	ENDIF
401	END DO
402	ind ( ji, jj ) = jk_in
403	e3_bak ( ji, jj ) = e3t_in ( ji, jj, jk_in ) ! store the value of the original layer thickness
404	e3t_in ( ji, jj, jk_in) = e3t_in ( ji, jj, jk_in ) - ( cff - hmax )
405	END DO
406	END DO
407	!
408	IF (Len_Trim(var_name) == 0) THEN
409	!/
410	! Interpolation of potential temperature TPOT
411	CALL zinterp ( jpi, jpj, jpka, jpka_in, ind, &
412	& tpot, e3t_in, e3_bak, e3t, varout, jp_weno )
413
414	varout(:,:,1,1) = varout(:,:,2,1)
415
416	IF (ln_c1d) THEN
417	DO jj=1,3
418	DO ji=1,3
419	varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
420	END DO
421	END DO
422	CALL Write_Ncdf_var ('tair', dimnames(1:4), out_file, varc1d, kt, 'float' )
423	ELSE
424	CALL Write_Ncdf_var ('tair', dimnames(1:4), out_file, varout, kt, 'float' )
425	ENDIF
426	!
427	! Interpolation of HUMI
428	CALL zinterp ( jpi, jpj, jpka, jpka_in, ind, &
429	& humi, e3t_in, e3_bak, e3t, varout, jp_weno )
430
431	!FL: dirty loop (possible issue in boundary conditions for WENO scheme)
432	DO jj = 1, jpj
433	DO ji = 1, jpi
434	DO jk = 2, jpka+1
435	varout(ji,jj,jk,1) = MAX(varout(ji,jj,jk,1),1.E-08) !<-- negative values in ECMWF
436	END DO
437	END DO
438	END DO
439
440	varout(:,:,1,1) = varout(:,:,2,1)
441
442	IF (ln_c1d) THEN
443	DO jj=1,3
444	DO ji=1,3
445	varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
446	END DO
447	END DO
448	CALL Write_Ncdf_var ('humi', dimnames(1:4), out_file, varc1d, kt, 'float' )
449	ELSE
450	CALL Write_Ncdf_var ('humi', dimnames(1:4), out_file, varout, kt, 'float' )
451	ENDIF
452
453	!
454	! Interpolate large-scale HPG or geostrophic wind
455	!
456	! Read HPG
457	IF (ln_hpg_frc) THEN
458	CALL Read_Ncdf_var ( varnames(8) , trim(cn_dir)//'/'//trim(file_in), tair, kt )
459	CALL Read_Ncdf_var ( varnames(9) , trim(cn_dir)//'/'//trim(file_in), humi, kt )
460	ENDIF
461	! Read geostrophic wind
462	IF (ln_geo_wnd) THEN
463	CALL Read_Ncdf_var ( varnames(10) , trim(cn_dir)//'/'//trim(file_in), tair, kt )
464	CALL Read_Ncdf_var ( varnames(11) , trim(cn_dir)//'/'//trim(file_in), humi, kt )
465	ENDIF
466	!
467	CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, tair( :,:,:,1 ) )
468	CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, humi( :,:,:,1 ) )
469	!
470	! Interpolation of geostrophic U
471	CALL zinterp ( jpi, jpj, jpka, jpka_in, ind, &
472	& tair, e3t_in, e3_bak, e3t, varout, jp_spln )
473	varout(:,:,1,1) = varout(:,:,2,1)
474
475	IF (ln_c1d) THEN
476	DO jj=1,3
477	DO ji=1,3
478	varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
479	END DO
480	END DO
481	CALL Write_Ncdf_var ( 'uhpg', dimnames(1:4), out_file, varc1d, kt, 'float' )
482	ELSE
483	IF (ln_hpg_frc) CALL Write_Ncdf_var ( varnames( 8), dimnames(1:4), out_file, varout, kt, 'float' )
484	IF (ln_geo_wnd) CALL Write_Ncdf_var ( varnames(10), dimnames(1:4), out_file, varout, kt, 'float' )
485	ENDIF
486	!
487	! Interpolation of geostrophic V
488	CALL zinterp ( jpi, jpj, jpka, jpka_in, ind, &
489	& humi, e3t_in, e3_bak, e3t, varout, jp_spln )
490	varout(:,:,1,1) = varout(:,:,2,1)
491
492	IF (ln_c1d) THEN
493	DO jj=1,3
494	DO ji=1,3
495	varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
496	END DO
497	END DO
498	CALL Write_Ncdf_var ( 'vhpg', dimnames(1:4), out_file, varc1d, kt, 'float' )
499	ELSE
500	IF (ln_hpg_frc) CALL Write_Ncdf_var ( varnames( 9), dimnames(1:4), out_file, varout, kt, 'float' )
501	IF (ln_geo_wnd) CALL Write_Ncdf_var ( varnames(11), dimnames(1:4), out_file, varout, kt, 'float' )
502	ENDIF
503	!
504	! Interpolation of total winds
505	!
506	! Read wind
507	CALL Read_Ncdf_var ( varnames(4) , trim(cn_dir)//'/'//trim(file_u), tair, kt )
508	CALL Read_Ncdf_var ( varnames(5) , trim(cn_dir)//'/'//trim(file_v), humi, kt )
509	CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, tair( :,:,:,1 ) )
510	CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, humi( :,:,:,1 ) )
511	!
512	! Interpolation of total U
513	CALL zinterp ( jpi, jpj, jpka, jpka_in, ind, &
514	& tair, e3t_in, e3_bak, e3t, varout, jp_spln )
515	varout(:,:,1,1) = varout(:,:,2,1)
516
517	IF(ln_c1d) THEN
518	DO jj=1,3
519	DO ji=1,3
520	varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
521	END DO
522	END DO
523	CALL Write_Ncdf_var ( 'uwnd', dimnames(1:4), out_file, varc1d, kt, 'float' )
524	ELSE
525	CALL Write_Ncdf_var ( 'uwnd', dimnames(1:4), out_file, varout, kt, 'float' )
526	ENDIF
527	!
528	! Interpolation of total V
529	CALL zinterp ( jpi, jpj, jpka, jpka_in, ind, &
530	& humi, e3t_in, e3_bak, e3t, varout, jp_spln )
531	varout(:,:,1,1) = varout(:,:,2,1)
532
533	IF(ln_c1d) THEN
534	DO jj=1,3
535	DO ji=1,3
536	varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
537	END DO
538	END DO
539	CALL Write_Ncdf_var ( 'vwnd', dimnames(1:4), out_file, varc1d, kt, 'float' )
540	ELSE
541	CALL Write_Ncdf_var ( 'vwnd', dimnames(1:4), out_file, varout, kt, 'float' )
542	ENDIF
543
544	ELSE ! var_name
545	!
546	ctrl = minval(pack([(ji,ji=1,size(outnames))],outnames==var_name))
547	var_file = filnames(ctrl)
548	!
549	!
550	! Interpolation of var_name
551	!
552	IF ( (var_name.NE."tair").AND.(var_name.NE."humi") ) THEN
553	! Read var_name
554	CALL Read_Ncdf_var ( varnames(ctrl) , trim(cn_dir)//'/'//trim(var_file), tair, kt )
555	CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, tair( :,:,:,1 ) )
556	! Interpolation of var_name
557	CALL zinterp ( jpi, jpj, jpka, jpka_in, ind, &
558	& tair, e3t_in, e3_bak, e3t, varout, jp_spln )
559	ELSE
560	! humi and tpot already read
561	IF (var_name.EQ."humi") tair = humi
562	IF (var_name.EQ."tair") tair = tpot
563	! Interpolation of var_name
564	CALL zinterp ( jpi, jpj, jpka, jpka_in, ind, &
565	& tair, e3t_in, e3_bak, e3t, varout, jp_weno )
566	END IF
567
568	varout(:,:,1,1) = varout(:,:,2,1)
569
570	IF(ln_c1d) THEN
571	DO jj=1,3
572	DO ji=1,3
573	varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
574	END DO
575	END DO
576	CALL Write_Ncdf_var ( outnames(ctrl), dimnames(1:4), out_file, varc1d, kt, 'float' )
577	ELSE
578	CALL Write_Ncdf_var ( outnames(ctrl), dimnames(1:4), out_file, varout, kt, 'float' )
579	ENDIF
580	!
581	ENDIF ! var_name
582	!
583	END DO ! kt
584	!
585	DEALLOCATE(zsurf,slp,tair,humi,varout)
586	IF (jpka_in.NE.nhym) DEALLOCATE(tmp_fullw,tmp_fullm)
587	!
588	STOP
589	!
590	END PROGRAM main

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

Download in other formats: