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_hdrown.F90 in utils/tools_ticket2362/ABL_TOOLS/src – NEMO

source: utils/tools_ticket2362/ABL_TOOLS/src/main_hdrown.F90 @ 13173

Last change on this file since 13173 was 13173, checked in by gsamson, 4 years ago

modify ABL_TOOLS directory structure to be maketools compliant; add ABL_TOOLS README.rst description; modify maketools to work with tools* directory names; ticket #2362
File size: 15.0 KB
Line 
1PROGRAM main
2   !!======================================================================
3   !!                     ***  PROGRAM main  ***
4   !!
5   !! ** Purpose : horizontal extrapolation ("drowning") of ECMWF dataset
6   !!
7   !!
8   !!======================================================================
9   !! History : 2016-10  (F. Lemarié)  Original code largely inspired by SOSIE (L. Brodeau & J.M. Molines)
10   !!   
11   !!----------------------------------------------------------------------
12   USE module_io       ! I/O routines
13   USE module_grid     ! compute input and output grids
14   !!
15   IMPLICIT NONE
16   !!----------------------------------------------------------------------
17   !!
18   !! 
19   !!
20   !!----------------------------------------------------------------------
21   !
22   INTEGER                                  :: ji,jj,jk,kt,kv,jjp1,jjm1,jip1,jim1
23   INTEGER                                  :: jpka,jpvar            ! number of vertical levels for input and target grids
24   INTEGER                                  :: jpi , jpj             ! number of grid points in x and y directions                 
25   INTEGER                                  :: jptime,ctrl,niter,status,ncid
26   INTEGER                                  :: ioerr
27   INTEGER, PARAMETER                       :: stdout   = 6
28   INTEGER, PARAMETER                       :: max_iter = 50
29   !!
30   REAL(8)                                  :: cff,cff1,cff2
31   !! 
32   REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: varin
33   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: tmask, tmask2, ff_t
34   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: mask_coast, maskv              ! land-sea mask   
35   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: varold,varnew   
36   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: tmp1d, lat
37   !!
38   CHARACTER(len=500)                       :: file_u,file_v,file_hpg,file_geos ! ERAi files containing wind components
39   CHARACTER(len=500)                       :: file_t,file_q, file_m   ! ERAi files containing tracers and mask
40   CHARACTER(len=500)                       :: file_z,file_p,cn_dir    ! ERAi files containing surface geopot and pressure
41   CHARACTER(len=500)                       :: abl_file, grd_file, drwn_file, out_file
42   CHARACTER(len=500)                       :: namelistf             
43   CHARACTER(len=500)                       :: argument   
44   CHARACTER(len= 20), DIMENSION(4)         :: dimnames
45   CHARACTER(:), ALLOCATABLE, DIMENSION(:)  :: varnames
46   CHARACTER(6)                             :: mask_var            ! name of mask variable in file_m file 
47   CHARACTER(6)                             :: var_name
48   !!
49   LOGICAL                                  :: ln_read_zsurf      ! read surface geopotential or not
50   LOGICAL                                  :: ln_read_mask       ! read land-sea mask or not
51   LOGICAL                                  :: ln_perio_latbc     ! use periodic BC along the domain latitudinal edges (for global data) or use zero-gradient BC (for regional data)
52   LOGICAL                                  :: ln_c1d             ! output only a single column in output file
53   LOGICAL                                  :: ln_hpg_frc         ! compute horizontal pressure gradient
54   LOGICAL                                  :: ln_geo_wnd         ! compute goestrophic wind components
55   LOGICAL                                  :: ln_slp_smth        ! apply gibbs oscillation filetring on mean sea level pressure
56   LOGICAL                                  :: ln_drw_smth        ! apply gibbs oscillation filetring on mean sea level pressure
57   LOGICAL                                  :: ln_slp_log         ! log(sea-level pressure) or sea-level pressure
58   LOGICAL                                  :: ln_lsm_land        ! if T mask is 1 over land and 0 over ocean if F it is the other way around
59   INTEGER                                  :: ptemp_method       ! way to compute potential temperature
60   !!
61   REAL(8), PARAMETER     :: omega =   7.292116e-05
62   REAL(8), PARAMETER     :: rad   =   3.141592653589793 / 180.
63
64
65   !!---------------------------------------------------------------------
66   !! List of variables read in the namelist file
67   NAMELIST/nml_out/    grd_file, abl_file, drwn_file, var_name
68   NAMELIST/nml_opt/    ptemp_method, ln_slp_log, ln_slp_smth, ln_read_mask, ln_perio_latbc, &
69      &                 ln_hpg_frc, ln_geo_wnd, ln_c1d, ln_read_zsurf, ln_lsm_land, ln_drw_smth
70   NAMELIST/nml_fld/    cn_dir, file_u, file_v, file_t,              &
71      &                 file_q, file_z, file_p, file_hpg, file_geos, &
72      &                 file_m, mask_var
73   !!
74   !! get the namelist file name
75   CALL get_command_argument(1, argument, ctrl, status)
76   !
77   SELECT CASE(status)
78   CASE(0)
79      namelistf = trim(argument)
80   CASE(-1)
81      WRITE(stdout,*) "### Error: file name too long"
82      STOP
83   CASE DEFAULT
84      namelistf = 'namelist_abl_tools'
85   END SELECT
86   !!---------------------------------------------------------------------
87
88
89   !!---------------------------------------------------------------------   
90   !! read namelist variables
91   ctrl = 0
92   OPEN(50, file=namelistf, status='old', form='formatted', access='sequential', iostat=ioerr)
93   IF (ioerr /= 0) ctrl = ctrl + 1   
94   READ(50,nml_opt, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
95   READ(50,nml_fld, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
96   READ(50,nml_out, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
97   
98   IF (ctrl > 0) then
99      WRITE(stdout,*) "### E R R O R while reading namelist file '",trim(namelistf),"'"
100      WRITE(stdout,*) " ctrl = ",ctrl
101      STOP
102   ELSE
103      WRITE(stdout,*) " Namelist file ",trim(namelistf)," OK "   
104   END IF
105   !!-------------------------------------------------------------------------------------   
106   !! list of variables to treat
107   !!
108   !! get the variable name
109   CALL get_command_argument( 2, argument, ctrl, status)
110   SELECT CASE(status)
111   CASE(0)
112      var_name = trim(argument)
113   CASE(-1)
114      WRITE(stdout,*) "### Error: file name too long"
115      STOP
116   END SELECT
117   WRITE(stdout,*) "var_name: ", trim(var_name), Len_Trim(var_name)
118   !
119   IF (Len_Trim(var_name) == 0) THEN
120     jpvar    = 7
121     allocate(character(len=20) :: varnames(jpvar) )
122     varnames(1:5) = [character(len=4) :: 'slp', 'uwnd', 'vwnd', 'tair', 'humi' ]
123     IF (ln_hpg_frc) varnames(6:7) = [character(len=4) :: 'uhpg', 'vhpg' ]
124     IF (ln_geo_wnd) varnames(6:7) = [character(len=4) :: 'ugeo', 'vgeo' ]
125     IF( .NOT. Var_Existence( varnames(5), trim(cn_dir)//'/'//trim(abl_file) ) ) jpvar = jpvar - 1
126     IF( .NOT. Var_Existence( varnames(6), trim(cn_dir)//'/'//trim(abl_file) ) ) jpvar = jpvar - 1 
127   ELSE
128     jpvar     = 1 
129     allocate(character(len=20) :: varnames(jpvar) )
130     varnames  = var_name
131     abl_file  = trim(var_name)//'_'//trim( abl_file)
132     drwn_file = trim(var_name)//'_'//trim(drwn_file)
133   END IF
134   WRITE(stdout,*) "Number of variables to treat : ",jpvar   
135   !!-------------------------------------------------------------------------------------   
136   !! read the dimensions for the input files
137   CALL Read_Ncdf_dim ( 'jpka' , trim(cn_dir)//'/'//trim(abl_file), jpka    ) 
138   CALL Read_Ncdf_dim ( 'time' , trim(cn_dir)//'/'//trim(abl_file), jptime  ) 
139   CALL Read_Ncdf_dim ( 'lon'  , trim(cn_dir)//'/'//trim(abl_file), jpi     )     
140   CALL Read_Ncdf_dim ( 'lat'  , trim(cn_dir)//'/'//trim(abl_file), jpj     )       
141   WRITE(stdout,*) "jpka, jptime, jpi, jpj: ", jpka, jptime, jpi, jpj
142   !
143   !!---------------------------------------------------------------------
144
145
146   !!---------------------------------------------------------------------
147   !! allocate arrays 
148   ALLOCATE(varin (1:jpi,1:jpj,1:jpka,1))
149   ALLOCATE(tmask (1:jpi,1:jpj    ))
150   ALLOCATE(tmask2(1:jpi,1:jpj    ))
151   ALLOCATE(maskv      (1:jpi,1:jpj )) 
152   ALLOCATE(mask_coast (1:jpi,1:jpj ))     
153   ALLOCATE(varold     (1:jpi,1:jpj )) 
154   ALLOCATE(varnew     (1:jpi,1:jpj ))     
155
156   !!---------------------------------------------------------------------
157   !! Read the mask and remove some closed seas
158   IF (ln_read_mask) THEN
159     CALL init_atm_mask(jpi,jpj,trim(cn_dir)//'/'//trim(file_m),trim(mask_var),ln_lsm_land, tmask )
160     CALL Read_Ncdf_var ( 'tmask' , trim(cn_dir)//'/'//trim(abl_file), tmask2(:,:) )
161   ELSE
162     tmask(:,:)  = 1.
163     tmask2(:,:) = 1.
164   END IF
165   !!   
166
167
168
169   !!---------------------------------------------------------------------   
170   !! create output file
171   !!
172   out_file = trim(cn_dir)//'/'//trim(drwn_file)
173   CALL Init_output_File ( jpi, jpj, jpka-1, trim(cn_dir)//'/'//trim(abl_file), out_file, tmask(:,:) )
174
175   !!---------------------------------------------------------------------   
176   !! Initialize the name of the dimensions for the result of the drowning
177   !!
178   dimnames(1) = 'lon'
179   dimnames(2) = 'lat'       
180   dimnames(3) = 'jpka'
181   dimnames(4) = 'time' 
182
183   CALL Write_Ncdf_var( 'tmask', dimnames(1:2), trim(out_file), tmask2, 'float' )
184
185   !!---------------------------------------------------------------------
186   ! Read time variable
187   ALLOCATE(tmp1d (1:jptime))   
188   ALLOCATE(lat   (1:jpj))   
189   CALL Read_Ncdf_var ( 'time', trim(cn_dir)//'/'//trim(file_t), tmp1d )   
190   CALL Read_Ncdf_var ( 'lat' , trim(cn_dir)//'/'//trim(file_t), lat )      !<-- latitude 
191   !!---------------------------------------------------------------------     
192
193   ! force drowning in the equatorial band with geo wind
194   IF (ln_geo_wnd) THEN
195     ALLOCATE( ff_t(1:jpi,1:jpj) )
196     DO jj = 1, jpj
197        DO ji = 1, jpi
198           ff_t(ji,jj) = 2. * omega * SIN( rad * lat(jj) )
199           IF (abs(ff_t(ji,jj)) < 2.5e-5) tmask2(ji,jj) = 0. 
200        END DO
201     END DO   
202   END IF
203
204   !===========
205   DO kv = 1,jpvar
206   !===========
207      DO kt = 1,jptime
208      !===========
209
210         IF( kv == 1 ) THEN
211           CALL Write_Ncdf_var( 'time', dimnames(4:4), trim(out_file), tmp1d(kt:kt), kt, 'double' )
212         ENDIF
213
214         ! Read variable to treat
215         print*,'Treat variable ',trim(varnames(kv)),' at time level ',kt
216         IF (trim(varnames(kv)).EQ."slp") THEN
217           CALL Read_Ncdf_var ( trim(varnames(kv)), trim(cn_dir)//'/'//trim(file_hpg), varin(:,:,1,1),  kt )
218         ELSE
219           CALL Read_Ncdf_var ( trim(varnames(kv)), trim(cn_dir)//'/'//trim(abl_file), varin,  kt )
220         END IF
221
222
223         !===========
224         DO jk = 1,jpka
225         !===========
226
227            IF ( (trim(varnames(kv))=="uwnd").OR.(trim(varnames(kv))=="vwnd").OR. &
228               & (trim(varnames(kv))=="ugeo").OR.(trim(varnames(kv))=="vgeo").OR. &
229               & (trim(varnames(kv))=="uhpg").OR.(trim(varnames(kv))=="vhpg") ) THEN
230                 maskv (1:jpi,1:jpj     ) = tmask2(1:jpi,1:jpj )
231            ELSE
232                 maskv (1:jpi,1:jpj     ) = tmask(1:jpi,1:jpj )
233            END IF
234
235
236            varold (1:jpi,1:jpj) = varin(1:jpi,1:jpj,jk,1)
237            varnew (1:jpi,1:jpj) = varin(1:jpi,1:jpj,jk,1)           
238            mask_coast(1:jpi,1:jpj) = 0.
239
240            !$$$$$$$$$$$$$$$$
241            DO niter = 1, max_iter
242            !$$$$$$$$$$$$$$$$                 
243
244               varold(1:jpi,1:jpj) = varnew (1:jpi,1:jpj)
245
246               IF ( .NOT. (ANY(maskv == 0))  ) THEN
247                  EXIT
248               END IF
249
250               !+++++++++++++++++++++++++++++++++++
251               ! Build mask_coast
252               DO jj = 1,jpj
253                  DO ji = 1,jpi
254
255                     IF (ln_perio_latbc) THEN
256                       jip1 = ji + 1; if(ji==jpi) jip1 = 1
257                       jim1 = ji - 1; if(ji==  1) jim1 = jpi
258                     ELSE
259                       jip1 = ji + 1; if(ji==jpi) jip1 = jpi-1
260                       jim1 = ji - 1; if(ji==  1) jim1 = 2
261                     END IF
262                     jjp1 = jj + 1; if(jj==jpj) jjp1 = jpj-1
263                     jjm1 = jj - 1; if(jj==  1) jjm1 = 2
264
265                     cff  = (1.-maskv(ji,jj))
266                     cff1 = maskv(ji,jjp1)+maskv(jip1,jj)+maskv(jim1,jj)+maskv(ji,jjm1) 
267                     mask_coast(ji,jj) = cff * cff1 
268                     IF( mask_coast(ji,jj) > 0. ) mask_coast(ji,jj) = 1.
269
270                  END DO
271               END DO   
272               !+++++++++++++++++++++++++++++++++++         
273
274               DO jj=1,jpj
275                  DO ji=1,jpi
276
277                     IF ( mask_coast(ji,jj) == 1. ) THEN
278
279                        IF (ln_perio_latbc) THEN
280                          jip1 = ji + 1; if(ji==jpi) jip1 = 1
281                          jim1 = ji - 1; if(ji==  1) jim1 = jpi
282                        ELSE
283                          jip1 = ji + 1; if(ji==jpi) jip1 = jpi-1
284                          jim1 = ji - 1; if(ji==  1) jim1 = 2
285                        END IF
286                        jjp1 = jj + 1; if(jj==jpj) jjp1 = jpj-1
287                        jjm1 = jj - 1; if(jj==  1) jjm1 = 2
288
289                        cff  = maskv(jim1,jjm1)+maskv(jim1,jj)+maskv(jim1,jjp1) &
290                           & + maskv(ji  ,jjm1)+               maskv(ji  ,jjp1) &
291                           & + maskv(jip1,jjm1)+maskv(jip1,jj)+maskv(jip1,jjp1)
292
293                        varnew(ji,jj) = (1./cff)*(                                   &
294                           &                    varold(jim1,jjm1)*maskv(jim1,jjm1)   &
295                           &              +     varold(jim1,jj  )*maskv(jim1,jj  )   &
296                           &              +     varold(jim1,jjp1)*maskv(jim1,jjp1)   &
297                           &              +     varold(ji  ,jjm1)*maskv(ji  ,jjm1)   &
298                           &              +     varold(ji  ,jjp1)*maskv(ji  ,jjp1)   &               
299                           &              +     varold(jip1,jjm1)*maskv(jip1,jjm1)   &
300                           &              +     varold(jip1,jj  )*maskv(jip1,jj  )   &
301                           &              +     varold(jip1,jjp1)*maskv(jip1,jjp1)   )
302
303                     END IF
304
305                  END DO
306               END DO   
307               !+++++++++++++++++++++++++++++++++++
308               maskv(1:jpi,1:jpj) = maskv(1:jpi,1:jpj) + mask_coast(1:jpi,1:jpj)
309               !+++++++++++++++++++++++++++++++++++                   
310           !$$$$$$$$$$$$$$$$
311           END DO
312           !$$$$$$$$$$$$$$$$
313
314           ! SMOOTHING AFTER DROWNING
315           IF (ln_drw_smth) THEN
316             maskv (1:jpi,1:jpj     ) = 1.
317           ELSE
318             !!only over continent
319             maskv (1:jpi,1:jpj     ) = 1. - tmask(1:jpi,1:jpj  )
320           END IF       
321
322           CALL smooth_field( jpi, jpj, varnew, maskv, 3 )
323           varin(1:jpi,1:jpj,jk,1) = varnew(1:jpi,1:jpj)
324
325           IF (trim(varnames(kv)).EQ."slp") EXIT
326
327        !===========
328        END DO ! jk
329        !===========
330
331        IF (trim(varnames(kv)).EQ."slp") THEN
332          CALL Write_Ncdf_var ( trim(varnames(kv)), dimnames(1:4), trim(cn_dir)//'/'//drwn_file, varin(:,:,1,1), kt, 'float' )
333        ELSE
334          CALL Write_Ncdf_var ( trim(varnames(kv)), dimnames(1:4), trim(cn_dir)//'/'//drwn_file, varin(:,:,:,:), kt, 'float' )
335        END IF
336
337      !===========
338      END DO ! time
339      !===========
340
341   !===========
342   END DO ! variable
343   !===========
344
345   STOP
346   !
347END PROGRAM main
Note: See TracBrowser for help on using the repository browser.