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

Last change on this file since 11589 was 11589, checked in by gsamson, 12 months 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 
1PROGRAM 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   !
590END PROGRAM main
Note: See TracBrowser for help on using the repository browser.