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

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

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

dev_r11265_ABL : see #2131

  • ABL_TOOLS first working version (README empty and arch files ignored for now)
File size: 34.6 KB
Line 
1MODULE module_grid
2   !!======================================================================
3   !!                   ***  MODULE  module_grid  ***
4   !! ABL utilities to define and store vertical grids
5   !!=====================================================================
6   !! History : 2016-10  (F. Lemarié)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   FUNCTIONS   : None   
11   !!   SUBROUTINES : get_atm_grid, init_atm_mask, get_pot_temp
12   !!                 Write_Grid_File, Init_output_File, init_target_grid
13   !!                 flip_vert_dim, smooth_field
14   !!----------------------------------------------------------------------
15   IMPLICIT NONE
16
17CONTAINS
18
19   SUBROUTINE get_atm_grid( jpi, jpj, jpka, slp, temp, humi, Aw, Bw, &   ! in
20      &                                                      e3t, ghw )   ! out
21      !!--------------------------------------------------------------------------
22      !!                    ***  ROUTINE get_atm_grid  ***
23      !!                   
24      !! ** Purpose : compute layer thickness and altitude of interfaces   
25      !!              of the ECMWF atmospheric grid
26      !!
27      !! ** Method  :
28      !!        (1) recompute the pressure levels thanks to the sea-level pressure
29      !!        (2) use the hydrostatic relation to convert pressure into altitudes
30      !!
31      !!---------------------------------------------------------------------------     
32      INTEGER, INTENT(in   ) ::   jpi, jpj
33      INTEGER, INTENT(in   ) ::   jpka
34      REAL(8), INTENT(in   ) ::   slp   ( 1:jpi, 1:jpj         )
35      REAL(8), INTENT(in   ) ::   temp  ( 1:jpi, 1:jpj, 1:jpka )     
36      REAL(8), INTENT(in   ) ::   humi  ( 1:jpi, 1:jpj, 1:jpka ) 
37      REAL(8), INTENT(in   ) ::   Aw    (               0:jpka )
38      REAL(8), INTENT(in   ) ::   Bw    (               0:jpka )                 
39      REAL(8), INTENT(  out) ::   e3t   ( 1:jpi, 1:jpj, 1:jpka ) 
40      REAL(8)                ::   ghw   ( 1:jpi, 1:jpj, 0:jpka )                           
41      REAL(8), PARAMETER     ::   g    =   9.80665
42      REAL(8), PARAMETER     ::   Rd   = 287.058   
43      REAL(8), PARAMETER     ::   zvir =   0.609133 
44      REAL(8), PARAMETER     ::   ig   =   1./g
45      !!
46      INTEGER                :: ji,jj,jk
47      REAL(8)                :: tv,ph(0:jpka)
48      !
49      DO jj = 1, jpj
50         DO ji = 1, jpi
51
52            DO jk=0,jpka
53               ph(jk) = Aw( jk ) + Bw( jk ) * slp( ji, jj )    !<-- Pa
54            END DO
55           
56            DO jk=jpka,1,-1
57              tv      = temp( ji, jj, jk ) * ( 1. + zvir*humi( ji, jj, jk ) )  !<-- Virtual temperature
58              e3t ( ji, jj, jk   ) =  ig*( Rd * tv * log( ph( jk ) / ph( jk-1 ) ) )
59              ghw ( ji, jj, jk-1 ) =  e3t( ji, jj, jk ) + ghw( ji, jj, jk )
60            END DO 
61
62         END DO
63      END DO
64      !
65   END SUBROUTINE get_atm_grid
66
67
68
69
70   SUBROUTINE init_atm_mask( jpi, jpj, mask_file, mask_name, ln_lsm_land, tmask )
71      USE module_io
72      !!---------------------------------------------------------------------
73      !!                    ***  ROUTINE INIT_atm_MASK  ***
74      !!                   
75      !! ** Purpose : extract the land/sea mask and remove isolated sea points
76      !!
77      !! ** Method  : mask is 1 over the ocean and 0 over land
78      !!
79      !!----------------------------------------------------------------------         
80      INTEGER, INTENT(in   ) ::   jpi, jpj 
81      LOGICAL, INTENT(in   ) ::   ln_lsm_land
82      REAL(8)                ::   tmask  ( 1:jpi, 1:jpj, 1 ) 
83      INTEGER                ::   jj, ji       
84      INTEGER                ::   status,ncid,varid 
85      CHARACTER(len = *  )   ::   mask_file, mask_name
86      REAL(8)                ::   cff
87       
88      ! Read land-sea mask variable
89      status    = nf90_open(trim(mask_file),NF90_NOWRITE,ncid)     
90      status    = nf90_inq_varid(ncid,mask_name,varid)       
91      status    = nf90_get_var(ncid,varid,tmask,start=(/1,1,1/))                                             
92      status    = nf90_close(ncid) 
93      ! invert the mask (1 over the ocean, 0 over land if ln_lsm_land)
94      IF(ln_lsm_land) THEN
95         DO jj=1,jpj
96            DO ji=1,jpi
97               IF( tmask(ji,jj,1) <= 0. ) THEN
98                  tmask(ji,jj,1) = 1.     !! Ocean points
99               ELSE
100                  tmask(ji,jj,1) = 0.     !! Land points     
101               END IF
102            END DO
103         END DO 
104      ENDIF
105      ! remove some closed seas
106      DO jj=2,jpj-1
107         DO ji=2,jpi-1
108            cff = MAX( tmask(ji+1,jj  ,1),tmask(ji-1,jj  ,1), &
109               &       tmask(ji  ,jj+1,1),tmask(ji  ,jj-1,1), &
110               &       tmask(ji+1,jj+1,1),tmask(ji-1,jj-1,1), &
111               &       tmask(ji+1,jj-1,1),tmask(ji-1,jj+1,1) )     
112            IF( tmask( ji, jj, 1 ) .gt. 0.5 .and. cff .lt. 0.5 ) THEN
113               tmask( ji, jj, 1 ) = 0.
114            END IF
115         END DO
116      END DO
117      WRITE(*,*)' init_atm_mask: ',mask_name,' in ',mask_file, ' OK'
118      !!----------------------------------------------------------------------
119      !
120   END SUBROUTINE init_atm_mask
121
122
123       
124   SUBROUTINE get_pot_temp( jpi, jpj, jpka, slp, temp, Aw, Bw, tmask, method, hum1, z1 )         
125      !!---------------------------------------------------------------------
126      !!                    ***  ROUTINE get_pot_temp  ***
127      !!                   
128      !! ** Purpose :  compute the potential temperature based on the
129      !!               absolute temperature and the sea level pressure
130      !!
131      !! ** Method  : five different ways are implemented depending on the
132      !!              value of 'method'
133      !!       (0) potential temperature = absolute temperature (not recommended)
134      !!       (1) potential temperature is computed using a local reference
135      !!                                pressure equal to the sea-level-pressure
136      !!       (2) potential temperature is computed only on a perturbation of
137      !!                                  the absolute temperature around t0
138      !!       (3) a local reference pressure is used consistently with AEROBULK gamma_moist
139      !!       (4) a constant global reference pressure is used
140      !!---------------------------------------------------------------------- 
141      INTEGER, INTENT(in   ) ::   jpi, jpj
142      INTEGER, INTENT(in   ) ::   jpka, method
143      REAL(8), INTENT(in   ) ::   slp   ( 1:jpi, 1:jpj         )
144      REAL(8), INTENT(in   ) ::   tmask ( 1:jpi, 1:jpj         )   
145      REAL(8), INTENT(in   ) ::   z1    ( 1:jpi, 1:jpj         )           
146      REAL(8), INTENT(inout) ::   temp  ( 1:jpi, 1:jpj, 1:jpka ) 
147      REAL(8), INTENT(in   ) ::   hum1  ( 1:jpi, 1:jpj         )       
148      REAL(8), INTENT(in   ) ::   Aw    (               0:jpka )
149      REAL(8), INTENT(in   ) ::   Bw    (               0:jpka )                         
150      REAL(8), PARAMETER     ::   grav   =   9.80665
151      REAL(8), PARAMETER     ::   R_dry  =   287.058 
152      REAL(8), PARAMETER     ::   R_vap  =   461.495
153      REAL(8), PARAMETER     ::   Cp_dry =   1005.             
154      REAL(8), PARAMETER     ::   cevap  =   2.5E+06       
155      REAL(8), PARAMETER     ::   reps0  =   R_dry / R_vap       
156      REAL(8), PARAMETER     ::   gamma  =   2./7.
157      REAL(8)                ::   pres0, gamma_moist, zrv, zirt     
158      REAL(8), PARAMETER     ::   t0     =      288.     !<-- K
159      !!
160      INTEGER                :: ji,jj,jk
161      REAL(8)                :: pres,ph(0:jpka), cff         
162      !!
163      SELECT CASE ( method ) 
164      CASE(0) 
165         RETURN
166      CASE(1)         
167         DO jj = 1, jpj
168            DO ji = 1, jpi
169
170               IF( tmask(ji,jj) .gt. 0.5 ) THEN
171                  DO jk=0,jpka
172                     ph(jk) = Aw( jk ) + Bw( jk ) * slp( ji, jj )
173                  END DO           
174
175                  pres0 = 0.5*(ph(jpka)+ph(jpka-1))
176           
177                  DO jk=1,jpka
178                     pres = 0.5*(ph(jk)+ph(jk-1))
179                     cff  = ( pres0/pres )**gamma
180                     temp( ji, jj, jk ) = cff * temp( ji, jj, jk )
181                  END DO
182               ELSE
183                  temp( ji, jj, 1:jpka ) = 273.15
184               END IF
185            END DO
186         END DO       
187      CASE(2)
188         pres0 =   100900.     
189         DO jj = 1, jpj
190            DO ji = 1, jpi
191               IF( tmask(ji,jj) .gt. 0.5 ) THEN
192                  DO jk=0,jpka
193                     ph(jk) = Aw( jk ) + Bw( jk ) * slp( ji, jj )
194                  END DO           
195           
196                  DO jk=1,jpka
197                     pres = 0.5*(ph(jk)+ph(jk-1))
198                     cff  = ( pres0/pres )**gamma
199                     temp( ji, jj, jk ) = cff * ( temp( ji, jj, jk ) - t0 ) + t0
200                  END DO
201               ELSE
202                  temp( ji, jj, 1:jpka ) = 273.15
203               END IF
204            END DO
205         END DO 
206      CASE(3)
207         DO jj = 1, jpj
208            DO ji = 1, jpi     
209               IF( tmask(ji,jj) .gt. 0.5 ) THEN
210                  DO jk=0,jpka
211                     ph(jk) = Aw( jk ) + Bw( jk ) * slp( ji, jj )
212                  END DO 
213                  !! compute gamma_moist consistently with AEROBULK
214                  zrv         = hum1( ji, jj ) / ( 1. - hum1( ji, jj ) )
215                  zirt        = 1. / ( R_dry * temp( ji, jj, jpka ) )
216                  gamma_moist = grav * ( 1. + cevap*zrv*ziRT )   & 
217                     &     / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/temp( ji, jj, jpka ) )       
218                  !! pressure at z = z1
219                  pres        = 0.5*(ph(jpka)+ph(jpka-1))                 
220                  pres0       = pres * ( 1. + gamma_moist * z1(ji,jj) / temp( ji, jj, jpka ) )**(1./gamma)
221                  DO jk=1,jpka
222                     pres = 0.5*(ph(jk)+ph(jk-1))
223                     cff  = ( pres0/pres )**gamma
224                     temp( ji, jj, jk ) = cff * temp( ji, jj, jk )
225                  END DO
226               ELSE
227                  temp( ji, jj, 1:jpka ) = 273.15
228               END IF
229            END DO
230         END DO           
231      CASE(4)
232         pres0 =   100900.     
233         DO jj = 1, jpj
234            DO ji = 1, jpi
235               IF( tmask(ji,jj) .gt. 0.5 ) THEN
236                  DO jk=0,jpka
237                     ph(jk) = Aw( jk ) + Bw( jk ) * slp( ji, jj )
238                  END DO           
239           
240                  DO jk=1,jpka
241                     pres = 0.5*(ph(jk)+ph(jk-1))
242                     cff  = ( pres0/pres )**gamma
243                     temp( ji, jj, jk ) = cff * temp( ji, jj, jk )
244                  END DO
245               ELSE
246                  temp( ji, jj, 1:jpka ) = 273.15
247               END IF
248            END DO
249         END DO       
250      END SELECT
251      !
252   END SUBROUTINE get_pot_temp
253
254
255   SUBROUTINE Write_Grid_File( jpka, ght, ghw, e3t, e3w, grd_file  ) 
256      !!
257      USE module_io
258      !!---------------------------------------------------------------------
259      !!                    ***  ROUTINE write_Grid_File  ***
260      !!                   
261      !! ** Purpose : write the ABL grid file 
262      !!
263      !! ** Method  : store the layer thickness and altitude of grid points
264      !!
265      !!----------------------------------------------------------------------
266      INTEGER,      INTENT(in   )  ::   jpka
267      REAL(8),      INTENT(in   )  ::   ght ( 1:jpka+1 )
268      REAL(8),      INTENT(in   )  ::   ghw ( 1:jpka+1 )     
269      REAL(8),      INTENT(in   )  ::   e3t ( 1:jpka+1 )
270      REAL(8),      INTENT(in   )  ::   e3w ( 1:jpka+1 )           
271      CHARACTER(*), INTENT(in   )  ::   grd_file
272      !!
273      INTEGER                      :: status, ncid
274      !!
275      status   = nf90_create( grd_file, NF90_WRITE, ncid )
276      status   = nf90_close ( ncid )   
277      !! 
278      Call Write_Ncdf_dim ( 'jpka' , grd_file, jpka+1 )       
279      Call Write_Ncdf_var ( 'ghw', 'jpka', grd_file, ghw, 'double' ) 
280      Call Write_Ncdf_var ( 'ght', 'jpka', grd_file, ght, 'double' )
281      Call Write_Ncdf_var ( 'e3t', 'jpka', grd_file, e3t, 'double' )
282      Call Write_Ncdf_var ( 'e3w', 'jpka', grd_file, e3w, 'double' )
283      !
284   END SUBROUTINE Write_Grid_File
285
286       
287
288
289
290   SUBROUTINE Init_output_File ( jpi, jpj, jpka, atm_file, abl_file, tmask  ) 
291      !!
292      USE module_io
293      !!---------------------------------------------------------------------
294      !!                    ***  ROUTINE Init_output_File  ***
295      !!                   
296      !! ** Purpose :  write longitude, latitude and mask in the output file
297      !!
298      !! ** Method  : define dimensions in the netcdf file
299      !!----------------------------------------------------------------------
300      INTEGER,                    INTENT(in   )  ::   jpi,jpj,jpka
301      CHARACTER(*),               INTENT(in   )  ::   atm_file     
302      CHARACTER(*),               INTENT(in   )  ::   abl_file
303      REAL(8)                                    ::   tmask(1:jpi,1:jpj)
304      !!
305      REAL(8), ALLOCATABLE, DIMENSION(:      )   ::   tmp1d
306      INTEGER                                    ::   status, ncid
307      CHARACTER(len= 20),DIMENSION(4)            :: dimnames
308      !!
309      status   = nf90_create( abl_file, NF90_WRITE, ncid )
310      status   = nf90_close ( ncid )   
311      !!     
312      CALL Write_Ncdf_dim ( 'lon'    , abl_file, jpi     )
313      CALL Write_Ncdf_dim ( 'lat'    , abl_file, jpj     )
314      CALL Write_Ncdf_dim ( 'jpka'   , abl_file, jpka+1  )         
315      CALL Write_Ncdf_dim ( 'time'   , abl_file, 0       )     
316      !
317      ALLOCATE( tmp1d( 1:jpi ) )
318      !
319      CALL Read_Ncdf_var  ( 'lon', atm_file, tmp1d )
320      CALL Write_Ncdf_var ( 'lon', 'lon', abl_file, tmp1d, 'double' ) 
321      !
322      DEALLOCATE( tmp1d )
323      !
324      ALLOCATE( tmp1d( 1:jpj ) )           
325      !
326      CALL Read_Ncdf_var ( 'lat', atm_file, tmp1d )   
327      CALL Write_Ncdf_var( 'lat', 'lat', abl_file, tmp1d, 'double' )     
328      !
329      DEALLOCATE( tmp1d )   
330      !
331      dimnames(1) = 'lon'
332      dimnames(2) = 'lat'
333      CALL Write_Ncdf_var( 'lsm', dimnames , abl_file, tmask, 'float' )
334      !
335   END SUBROUTINE Init_output_File
336       
337       
338   SUBROUTINE Init_output_File_c1d ( jpi, jpj, jpka, atm_file, abl_file, tmask, iloc, jloc )
339      !!
340      USE module_io
341      !!---------------------------------------------------------------------
342      !!                    ***  ROUTINE Init_output_File  ***
343      !!
344      !! ** Purpose :  write longitude, latitude and mask in the output file
345      !!
346      !! ** Method  : define dimensions in the netcdf file
347      !!----------------------------------------------------------------------
348      INTEGER,                    INTENT(in   )  ::   iloc,jloc,jpi,jpj,jpka
349      CHARACTER(*),               INTENT(in   )  ::   atm_file
350      CHARACTER(*),               INTENT(in   )  ::   abl_file
351      REAL(8)                                    ::   tmask(1:jpi,1:jpj)
352      !!
353      REAL(8), ALLOCATABLE, DIMENSION(:      )   ::   tmp1d
354      REAL(8)                                    ::   tmp1d_loc(3),tmp2d_loc(3,3)
355      INTEGER                                    ::   status, ncid
356      CHARACTER(len= 20),DIMENSION(4)            :: dimnames
357      !!
358      status   = nf90_create( abl_file, NF90_WRITE, ncid )
359      status   = nf90_close ( ncid )
360      !!
361      CALL Write_Ncdf_dim ( 'lon'    , abl_file, 3       )
362      CALL Write_Ncdf_dim ( 'lat'    , abl_file, 3       )
363      CALL Write_Ncdf_dim ( 'jpka'   , abl_file, jpka+1  )
364      CALL Write_Ncdf_dim ( 'time'   , abl_file, 0       )
365      !
366      ALLOCATE( tmp1d    ( 1:jpi ) )
367      !
368      CALL Read_Ncdf_var  ( 'lon', atm_file, tmp1d )
369      tmp1d_loc( 1 ) = tmp1d(iloc)
370      tmp1d_loc( 2 ) = tmp1d(iloc)
371      tmp1d_loc( 3 ) = tmp1d(iloc)
372      CALL Write_Ncdf_var ( 'lon', 'lon', abl_file, tmp1d_loc, 'double' )
373      print*,'1D column located at : '
374      print*,'longitude = ',tmp1d(iloc),' degree_east'
375      !
376      DEALLOCATE( tmp1d )
377      !
378      ALLOCATE( tmp1d( 1:jpj ) )
379      !
380      CALL Read_Ncdf_var ( 'lat', atm_file, tmp1d )
381      tmp1d_loc( 1 ) = tmp1d(jloc)
382      tmp1d_loc( 2 ) = tmp1d(jloc)
383      tmp1d_loc( 3 ) = tmp1d(jloc)
384      CALL Write_Ncdf_var( 'lat', 'lat', abl_file, tmp1d_loc, 'double' )
385      print*,'latitude = ',tmp1d(jloc),' degree_north'
386      !
387      DEALLOCATE( tmp1d )
388      !
389      tmp2d_loc(:,:) = tmask(iloc,jloc)       
390      dimnames(1) = 'lon'
391      dimnames(2) = 'lat'
392      CALL Write_Ncdf_var( 'lsm', dimnames , abl_file, tmp2d_loc, 'float' )
393      !
394   END SUBROUTINE Init_output_File_c1d       
395       
396
397   SUBROUTINE init_target_grid( jpka, ght, ghw, e3t, e3w, hmax, hc, theta_s,  &
398      &                                                        force_z1, z1 )
399      !!---------------------------------------------------------------------
400      !!                    ***  ROUTINE init_target_grid  ***
401      !!                   
402      !! ** Purpose : compute the layer thickness and altitude of grid points
403      !!              for the ABL model based on the namelist parameter values 
404      !!
405      !! ** Method  : depending on the logical 'force_z1' two methods are used
406      !!              (1) true -> the user chooses the value of the first vertical
407      !!                  grid point. A few Newton iterations are used to correct 
408      !!                  the value of the parameter theta_s to satisfy this constraint
409      !!              (2) false -> use the parameter values in the namelist to
410      !!                  compute the vertical grid
411      !!
412      !!----------------------------------------------------------------------
413      INTEGER, intent(in   )    :: jpka
414      REAL(8), intent(in   )    :: hmax,hc
415      REAL(8), intent(inout)    :: theta_s     
416      LOGICAL, intent(in   )    :: force_z1
417      REAL(8), intent(in   )    :: z1
418      !!
419      REAL(8), intent(  out)    :: ghw( 1:jpka+1 )
420      REAL(8), intent(  out)    :: ght( 1:jpka+1 )     
421      REAL(8), intent(  out)    :: e3w( 1:jpka+1 )
422      REAL(8), intent(  out)    :: e3t( 1:jpka+1 ) 
423      !!
424      REAL(8)                   ::  ds,cff,sc_w,sc_r,alpha,x
425      REAL(8)                   ::  fx,fxp
426      INTEGER                   ::  jk,maxiter,jiter
427      REAL(8), PARAMETER        :: tol = 1.E-12
428      !!
429      IF(force_z1) THEN
430        IF(z1.LT.10.) THEN
431          WRITE(*,*) " ERROR: z1 < 1st ECMWF level height (~10m)"
432          STOP
433        ELSE
434          !! Newton iterations to find the appropriate value of theta_s
435          maxiter = 1000
436          x       = theta_s
437          sc_r    = (float(1)-0.5)/float(jpka)
438          alpha   = (z1 - hc*sc_r) / (hmax - hc)
439          !
440          DO jiter=1,maxiter
441             fx   = (sinh(sc_r*x)/sinh(x))-alpha
442             fxp  = (sc_r*cosh(sc_r*x)-sinh(sc_r*x)*cosh(x)/sinh(x))/sinh(x)
443             IF( abs(fx) .lt. tol ) THEN
444                   exit
445             ENDIF
446             cff  = fx / fxp
447             x    = x - cff     
448          ENDDO
449          !
450          theta_s = x
451        END IF
452        !   
453      ENDIF
454      !
455      ds =1./float(jpka)
456      cff=(hmax-hc)/sinh(theta_s)
457      !
458      DO jk = jpka,1,-1
459        sc_w      = ds*float(jk)
460        ghw(jk+1) = hc*sc_w + cff*sinh(theta_s*sc_w)
461        sc_r      = ds*(float(jk)-0.5)
462        ght(jk+1) = hc*sc_r + cff*sinh(theta_s*sc_r)
463      END DO
464      !
465      ghw(1) = 0.
466      e3t(1) = 0.
467      ght(1) = 0.
468      !
469      DO jk = 2,jpka+1
470        e3t(jk) = ghw(jk)-ghw(jk-1)
471      END DO
472      !
473      DO jk=1,jpka
474        e3w(jk)  = ght(jk+1)-ght(jk)
475      END DO
476      !
477      e3w(jpka+1) = ghw(jpka+1) - ght(jpka+1)         
478      !
479      IF(force_z1) THEN   !++ print the new parameter values
480         print*,'*** Updated grid parameters'
481         print*,'theta_s = ',theta_s
482         print*,'hc      = ',hc     
483         print*,'hmax    = ',hmax         
484         print*,'ght(2)  = ',ght(2)
485      ENDIF     
486      !
487   END SUBROUTINE init_target_grid
488   
489   
490   
491   
492   
493
494   SUBROUTINE flip_vert_dim( kstr, kend, jpi, jpj, tabin )
495      !!---------------------------------------------------------------------
496      !!                    ***  ROUTINE flip_vert_dim  ***
497      !!                   
498      !! ** Purpose : flip the vertical axis of the array tabin so that   
499      !!              the vertical grid goes from k=kstr at the bottom
500      !!              of the ABL to k=kend at the top
501      !!----------------------------------------------------------------------
502      INTEGER, intent(in   )   :: kstr,kend
503      INTEGER, intent(in   )   :: jpi, jpj 
504      REAL(8), intent(inout)   :: tabin( 1:jpi, 1:jpj, kstr:kend )
505      !!
506      INTEGER                  :: ji,jj,jk,ks
507      REAL(8)                  :: tabfl(kstr:kend)
508      !
509      DO jj = 1,jpj
510         DO ji = 1,jpi
511            DO jk=kstr,kend
512               ks=(kend-jk)+kstr
513               tabfl(ks) = tabin( ji, jj, jk ) 
514            END DO
515            tabin( ji, jj, kstr:kend ) = tabfl( kstr:kend )
516         END DO
517      END DO 
518      !
519   END SUBROUTINE flip_vert_dim
520   
521 
522
523
524
525
526
527   SUBROUTINE smooth_field( jpi, jpj, varin, tmask, niter )         
528      !!---------------------------------------------------------------------
529      !!                    ***  ROUTINE smooth_field  ***
530      !!                   
531      !! ** Purpose : smooth the sea level pressure over the ocean
532      !!              to attenuate Gibbs oscillation
533      !! ** Method  : a 9-point isotropic laplacian filter is applied
534      !!              iteratively on ocean grid points only
535      !!
536      !! Proper treatment of the periodicity is still not yet implemented
537      !!
538      !!----------------------------------------------------------------------     
539      INTEGER, INTENT(in   ) ::   jpi, jpj, niter 
540      REAL(8), INTENT(in   ) ::   tmask  ( 1:jpi, 1:jpj )       
541      REAL(8), INTENT(inout) ::   varin  ( 1:jpi, 1:jpj )
542      INTEGER                :: ji,jj,nit
543      REAL(8)                :: smth_a,smth_b,umask,vmask 
544      REAL(8)                :: FX ( 0:jpi  , 0:jpj+1 )
545      REAL(8)                :: FE1( 0:jpi+1, 0:jpj   )
546      REAL(8)                :: FE ( 1:jpi  , 0:jpj   )       
547      !! 
548      !!=========================================================
549      !!   
550      ! Hanning filter
551      !smth_a = 1./8.
552      !smth_b = 1./4.     
553      ! 9-point isotropic laplacian filter
554      smth_a = 1./12.
555      smth_b = 3./16.     
556
557      FX ( 0:jpi  , 0:jpj+1 ) = 0.
558      FE1( 0:jpi+1, 0:jpj   ) = 0.
559      FE ( 1:jpi  , 0:jpj   ) = 0.
560
561      !!+++++++++
562      DO nit = 1,niter
563      !!+++++++++
564         DO jj=1,jpj             
565            DO ji=1,jpi-1         
566            umask = tmask(ji,jj)*tmask(ji+1,jj) 
567            FX ( ji, jj ) = ( varin( ji+1,jj ) - varin( ji  ,jj ) ) * umask 
568            END DO                   
569         END DO
570         FX(   0  , 1:jpj ) = FX( jpi-1, 1:jpj  )
571         FX  ( jpi, 1:jpj ) = FX(     1, 1:jpj  )     
572         FX( 0:jpi, 0     ) = FX( 0:jpi,     1  )
573         FX( 0:jpi, jpj+1 ) = FX( 0:jpi,     jpj)           
574         !!
575         DO jj=1,jpj-1                 
576            DO ji=1,jpi
577               vmask = tmask(ji,jj)*tmask(ji,jj+1) 
578               FE1( ji, jj ) = ( varin( ji, jj+1 ) - varin( ji  ,jj ) ) * vmask 
579           END DO
580         END DO
581         !!
582         FE1( 0      , 1:jpj-1 ) = FE1( jpi , 1:jpj-1 )
583         FE1( jpi+1  , 1:jpj-1 ) = FE1(   1 , 1:jpj-1 )     
584         FE1( 0:jpi+1,    0    ) = 0.
585         FE1( 0:jpi+1,    jpj  ) = 0. 
586         !!
587         DO jj=0,jpj                 
588            DO ji=1,jpi
589               FE ( ji, jj ) =   FE1( ji, jj   )                   &
590               & +  smth_a*(  FX ( ji, jj+1 ) + FX( ji-1, jj   )   &
591               &           -  FX ( ji, jj   ) - FX( ji-1, jj+1 ) )
592            END DO
593         END DO
594         !!
595         DO jj = 1, jpj
596            DO ji = 0,jpi
597               FX( ji, jj ) = FX( ji, jj )         &
598              &    + smth_a*( FE1( ji+1, jj ) + FE1( ji  , jj-1 )    &
599              &              -FE1( ji  , jj ) - FE1( ji+1, jj-1 ) )
600            END DO
601            DO ji = 1,jpi
602               varin( ji  ,jj ) = varin( ji  ,jj ) &
603     &         + tmask(ji,jj) * smth_b * (         &
604     &                  FX( ji, jj ) - FX( ji-1, jj )          &
605     &                 +FE( ji, jj ) - FE( ji, jj-1 )  )
606            END DO
607         END DO     
608         !!
609      !!+++++++++     
610      END DO
611      !!+++++++++   
612           
613   END SUBROUTINE smooth_field
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630   
631
632 
633
634
635
636
637!   SUBROUTINE DTV_Filter( jpi, jpj, varin, tmask, niter, time )         
638!      USE module_io
639!      !!---------------------------------------------------------------------
640!      !!                    ***  ROUTINE DTV_Filter  ***
641!      !!                   
642!      !! ** Purpose :   
643!      !!
644!      !! ** Method  :
645!      !!
646!      !! ** Action  :
647!      !!----------------------------------------------------------------------   
648!      INTEGER, INTENT(in   ) ::   jpi, jpj, niter, time
649!      REAL(8), INTENT(in   ) ::   tmask  ( 1:jpi, 1:jpj )       
650!      REAL(8), INTENT(inout) ::   varin  ( 1:jpi, 1:jpj )
651!      INTEGER                ::   jn, jj, ji, nit, nt_n, nt_a, nl2     
652!      REAL(8)                ::   var0,mean,sigma2,cff,lambda
653!      REAL(8)                ::   FX(0:jpi+1,0:jpj+1)
654!      REAL(8)                ::   FE(0:jpi+1,0:jpj+1)
655!      REAL(8)                ::   FL(0:jpi+1,0:jpj+1)
656!      REAL(8)                ::   FR(0:jpi+1,0:jpj+1) 
657!      REAL(8)                ::   wrk ( 1:jpi, 1:jpj, 2 )
658!      REAL(8)                ::   div ( 0:jpi+1, 0:jpj+1 ), diag               
659!      REAL(8), PARAMETER     ::   rsmall = 1.E-08
660!      REAL(8), PARAMETER     ::   rbig   = 1.E+14
661!      REAL(8)                ::   umask, vmask, fmask, wght(8), dt
662!      REAL(8)                ::   L2norm_n, L2norm_a
663!      REAL(8), ALLOCATABLE, DIMENSION(:      )   ::   tmp1d     
664!      INTEGER                :: status,ncid,varid1,varid2,dimid1,dimid2 
665!      !!=========================================================
666!      CHARACTER(len=  3  )   :: nn
667!      CHARACTER(len = 500)   :: erai_file, smth_file
668!      LOGICAL                :: ln_diag_smoothing
669!      ln_diag_smoothing = .true.
670!      IF( ln_diag_smoothing ) THEN
671!         erai_file = '/Users/florianlemarie/Documents/INRIA/SIMBAD/ERA_INTERIM_DECEMBRE_2015/phi_ml1_6h_erai_201512.nc'
672!         WRITE(nn,'(I3)') time
673!         DO ji=1,3
674!            IF(nn(ji:ji)==' ') nn(ji:ji)='0'
675!         END DO
676!     
677!         smth_file = 'smoothing_results_DTV_'//nn//'.nc' 
678!         status    = nf90_create( trim(smth_file) , NF90_WRITE, ncid )
679!         status    = nf90_close ( ncid )   
680!         !   
681!         CALL Write_Ncdf_dim ( 'lon'    , trim(smth_file) , jpi     )
682!         CALL Write_Ncdf_dim ( 'lat'    , trim(smth_file) , jpj     )
683!         !
684!         ALLOCATE( tmp1d( 1:jpi ) )
685!         CALL Read_Ncdf_var  ( 'lon', erai_file, tmp1d )
686!         CALL Write_Ncdf_var ( 'lon', 'lon', trim(smth_file), tmp1d, 'double' )
687!         DEALLOCATE( tmp1d )
688!         !
689!         ALLOCATE( tmp1d( 1:jpj ) )           
690!         CALL Read_Ncdf_var ( 'lat', erai_file, tmp1d )   
691!         CALL Write_Ncdf_var( 'lat', 'lat', trim(smth_file), tmp1d, 'double' )     
692!         DEALLOCATE( tmp1d )
693!         !         
694!         status = nf90_open(trim(smth_file),NF90_WRITE,ncid)     
695!         status = nf90_inq_dimid(ncid, 'lon', dimid1)
696!         status = nf90_inq_dimid(ncid, 'lat', dimid2)
697!         status = nf90_redef(ncid)     
698!         status = nf90_def_var(ncid,'varinp',nf90_double,(/dimid1,dimid2/),varid1)       
699!         status = nf90_def_var(ncid,'varout',nf90_double,(/dimid1,dimid2/),varid2)   
700!         status = nf90_enddef(ncid)
701!         !
702!         !
703!         status = nf90_put_var(ncid,varid1,varin)
704!      END IF
705!      !!=========================================================
706!      nit    = 0
707!      nl2    = 0
708!      L2norm_a = 0.
709!      L2norm_n = rbig     
710!      nt_n   = 1 + MOD( nit  , 2 )
711!      nt_a   = 1 + MOD( nit+1, 2 )     
712!      var0   = varin(  nint(0.5*jpi) , nint(0.5*jpj) ) 
713!      varin( 1:jpi, 1:jpj       ) = varin( 1:jpi, 1:jpj ) - var0
714!      wrk  ( 1:jpi, 1:jpj, nt_n ) = varin( 1:jpi, 1:jpj )
715!      !! Compute the mean first
716!      mean   = Get_Mean ( jpi , jpj,  varin, tmask       )
717!      sigma2 = Get_Vari ( jpi , jpj,  varin, tmask, mean )
718!      !!
719!      print*,'mean value over the ocean = ',mean + var0
720!      print*,'variance of the input field = ',sigma2     
721!      !!
722!      lambda = (1. / sigma2      )
723!      dt     =  1. / (16.+lambda)
724!      diag   =  0.
725!     
726!      !!>>>>>>>>>>>>>>>>
727!      DO nit = 1,niter
728!      !!>>>>>>>>>>>>>>>>
729!      FX(0:jpi+1,0:jpj+1) = 0.
730!      FE(0:jpi+1,0:jpj+1) = 0.
731!      FL(0:jpi+1,0:jpj+1) = 0.
732!      FR(0:jpi+1,0:jpj+1) = 0.
733!      !!
734!      DO jj = 1,jpj
735!         DO ji = 1,jpi-1
736!            umask        = tmask(ji,jj)*tmask(ji+1,jj) 
737!            FX( ji, jj ) = umask * ( wrk( ji+1, jj, nt_n ) - wrk( ji, jj, nt_n ) )
738!         END DO
739!      END DO 
740!      !!
741!      DO jj = 1,jpj-1
742!         DO ji = 1,jpi
743!            vmask        = tmask(ji,jj)*tmask(ji,jj+1)
744!            FE( ji, jj ) = vmask * ( wrk( ji, jj+1, nt_n ) - wrk( ji, jj, nt_n ) )
745!         END DO
746!      END DO       
747!      !!
748!      DO jj = 1,jpj-1
749!         DO ji = 1,jpi-1
750!            fmask        = tmask(ji,jj)*tmask(ji+1,jj+1)*diag
751!            FL( ji, jj ) = fmask * ( wrk( ji+1, jj+1, nt_n ) - wrk( ji  , jj, nt_n   ) )
752!            fmask        = tmask(ji+1,jj)*tmask(ji,jj+1)*diag       
753!            FR( ji, jj ) = fmask * ( wrk( ji  , jj+1, nt_n ) - wrk( ji+1, jj, nt_n ) )
754!         END DO
755!      END DO
756!      !!
757!      div( 0:jpi+1, 0:jpj+1 ) = rsmall
758!      !!
759!      DO jj = 1,jpj+1
760!         DO ji = 1,jpi+1           
761!            div(ji,jj) = MAX( sqrt(  FX(ji,jj)**2 + FX(ji-1,jj  )**2 + FE(ji,jj  )**2 + FX(ji,jj-1)**2          &
762!            &                      + FL(ji,jj)**2 + FL(ji-1,jj-1)**2 + FR(ji,jj-1)**2 + FR(ji-1,jj)**2 ), rsmall )
763!         END DO
764!      END DO
765!     
766!     
767!       
768!     
769!      DO jj=1,jpj
770!         DO ji=1,jpi
771!            IF( div(ji,jj) .eq. rsmall .or. tmask(ji,jj).lt.0.5 ) THEN
772!               wrk(ji,jj,nt_a) = wrk(ji,jj,nt_n)
773!            ELSE
774!               wght(1)     = 1. +  div(ji,jj) / div(ji+1,jj  )
775!               wght(2)     = 1. +  div(ji,jj) / div(ji-1,jj  )
776!               wght(3)     = 1. +  div(ji,jj) / div(ji  ,jj+1)                       
777!               wght(4)     = 1. +  div(ji,jj) / div(ji  ,jj-1)           
778!               wght(5)     = 1. +  div(ji,jj) / div(ji+1,jj+1)   
779!               wght(6)     = 1. +  div(ji,jj) / div(ji-1,jj+1)
780!               wght(7)     = 1. +  div(ji,jj) / div(ji-1,jj-1)                       
781!               wght(8)     = 1. +  div(ji,jj) / div(ji+1,jj-1)             
782!             
783!            wrk(ji,jj,nt_a) = wrk(ji,jj,nt_n) +  dt*(                   &
784!               &            + FX(ji  ,jj  ) * wght(1)                   &
785!               &            - FX(ji-1,jj  ) * wght(2)                   &
786!               &            + FE(ji  ,jj  ) * wght(3)                   &
787!               &            - FE(ji  ,jj-1) * wght(4)                   &
788!               &            + FL(ji  ,jj  ) * wght(5)                   &
789!               &            + FR(ji-1,jj  ) * wght(6)                   &
790!               &            - FL(ji-1,jj-1) * wght(7)                   &
791!               &            - FR(ji  ,jj-1) * wght(8)                   &
792!               &            - lambda * div(ji,jj) * ( wrk(ji,jj,nt_n) - varin(ji,jj) ) ) 
793!               
794!              IF( isnan(wrk(ji,jj,nt_a)) ) THEN
795!                  print*,'Nan in smoothing at iteration ',nit
796!                  print*,'at grid point ',ji,jj
797!              END IF   
798!              IF( abs(wrk(ji,jj,nt_a)) .gt. rbig ) THEN
799!                  print*,'Inf in smoothing at iteration ',nit
800!                  print*,'at grid point ',ji,jj                                   
801!              END IF
802!                 L2norm_a = L2norm_a + ( wrk(ji,jj,nt_n)-wrk(ji,jj,nt_a) )**2
803!              END IF                 
804!         END DO
805!      END DO
806!
807!
808!      mean   = Get_Mean ( jpi , jpj,  wrk(:,:,nt_a), tmask       )
809!      sigma2 = Get_Vari ( jpi , jpj,  wrk(:,:,nt_a), tmask, mean )
810!      !!
811!      print*,'mean value over the ocean = ',mean + var0
812!      print*,'variance of the input field = ',sigma2   
813!
814!     
815!      IF( L2norm_a .gt. L2norm_n ) THEN
816!         print*,'convergence after ',nit,' iterations'
817!!         EXIT
818!      END IF     
819!      L2norm_n = L2norm_a
820!      L2norm_a = 0.     
821!      nt_n = 1 + MOD( nit  , 2 )
822!      nt_a = 1 + MOD( nit+1, 2 )
823!           
824!      !!>>>>>>>>>>>>>>>>
825!      END DO
826!      !!>>>>>>>>>>>>>>>>
827!
828!       
829!      DO jj=1,jpj
830!         DO ji=1,jpi
831!            varin(ji,jj) = ( wrk(ji,jj,nt_n) + var0 )*tmask(ji,jj)
832!         END DO
833!      END DO
834!
835!      IF( ln_diag_smoothing ) THEN
836!         status = nf90_put_var(ncid,varid2,varin)
837!         status = nf90_close(ncid)
838!      END IF   
839!      !
840!   END SUBROUTINE dtv_filter
841!
842!
843!
844!   REAL(8) FUNCTION Get_Mean ( jpi , jpj, tabvar, tmask )
845!      !!---------------------------------------------------------------------
846!      !!                    ***  FUNCTION Get_Mean  ***
847!      !!                   
848!      !! ** Purpose : get the mean of the input field
849!      !!
850!      !!----------------------------------------------------------------------
851!      INTEGER,   INTENT(in)      :: jpi,jpj
852!      REAL(8),   INTENT(in)      :: tabvar(jpi,jpj)
853!      REAL(8),   INTENT(in)      :: tmask (jpi,jpj)
854!      INTEGER                    :: ji,jj,jn   
855!     
856!      !! Compute the mean first
857!      jn       = 0
858!      Get_Mean = 0.
859!      DO jj = 1, jpj
860!         DO ji = 1, jpi
861!            IF(tmask(ji,jj) .gt. 0.5) THEN
862!               Get_Mean = Get_Mean + tabvar( ji,jj ) 
863!               jn=jn+1
864!            END IF 
865!         END DO
866!      END DO     
867!      Get_Mean   = (1./jn)   * Get_Mean 
868!      !
869!   END FUNCTION Get_Mean
870!     
871!
872!
873!
874!
875!   REAL(8) FUNCTION Get_Vari ( jpi , jpj, tabvar, tmask, mean )
876!      !!---------------------------------------------------------------------
877!      !!                    ***  FUNCTION Get_Mean  ***
878!      !!                   
879!      !! ** Purpose : get the mean of the input field
880!      !!
881!      !!----------------------------------------------------------------------
882!      INTEGER,   INTENT(in)      :: jpi,jpj
883!      REAL(8),   INTENT(in)      :: tabvar(jpi,jpj)
884!      REAL(8),   INTENT(in)      :: tmask (jpi,jpj), mean
885!      INTEGER                    :: ji,jj,jn   
886!     
887!      jn       = 0
888!      Get_Vari = 0
889!      !!
890!      DO jj = 1, jpj
891!         DO ji = 1, jpi
892!            IF(tmask(ji,jj) .gt. 0.5) THEN
893!              Get_Vari = Get_Vari + ( tabvar(ji,jj) - mean )**2
894!              jn = jn+1
895!            ENDIF 
896!         END DO
897!      END DO
898!      Get_Vari   = (1./jn)   * Get_Vari
899!      !
900!   END FUNCTION Get_Vari
901
902
903   
904END MODULE module_grid
Note: See TracBrowser for help on using the repository browser.