source: utils/tools/ABL_TOOLS/main_uvg_hpg.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.7 KB
Line 
1PROGRAM 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   !
647END PROGRAM main
Note: See TracBrowser for help on using the repository browser.