3dc55d7c711eccc5329848f6c822f4188a9a5020
[MNH-git_open_source-lfs.git] / src / SURFEX / pgd_grid.F90
1 !SURFEX_LIC Copyright 1994-2014 Meteo-France 
2 !SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C  licence
3 !SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SURFEX_LIC for details. version 1.
5 !     #########
6       SUBROUTINE PGD_GRID(HPROGRAM,HFILE,HFILETYPE,OGRID,HGRID,KGRID_PAR,PGRID_PAR)
7 !     ##########################################################
8 !!
9 !!    PURPOSE
10 !!    -------
11 !!   Reads in namelist the grid type and parameters.
12 !!
13 !!    METHOD
14 !!    ------
15 !!   
16 !!    EXTERNAL
17 !!    --------
18 !!
19 !!
20 !!    IMPLICIT ARGUMENTS
21 !!    ------------------
22 !!
23 !!
24 !!    REFERENCE
25 !!    ---------
26 !!
27 !!    AUTHOR
28 !!    ------
29 !!
30 !!    V. Masson                   Meteo-France
31 !!
32 !!    MODIFICATION
33 !!    ------------
34 !!
35 !!    Original     01/2004
36 !!    E. Martin    10/2007 IGN grid
37 !!    M. Moge      05/02/2015 parallelization (using local sizes, GET_MEAN_OF_COORD_SQRT_ll, SET_NAM_GRID_CONF_PROJ_LOCAL) + MPPDB_CHECK
38 !!    M. Moge      01/03/2015 call SPLIT_GRID if CPROGRAM == 'PGD   ' + remove SET_NAM_GRID_CONF_PROJ_LOCAL
39 !!    M. Moge      01/03/2015 change in the input arguments of PGD_GRID_IO_INIT : passing IDXRATIO, IDYRATIO
40 !----------------------------------------------------------------------------
41 !
42 !*    0.     DECLARATION
43 !            -----------
44 !
45 USE MODD_SURFEX_MPI,     ONLY : NSIZE, NINDEX, NPIO, NRANK
46 USE MODD_SURFEX_OMP,     ONLY : NINDX2, NWORK, XWORK, XWORK2, XWORK3, &
47                                 NWORK_FULL, XWORK_FULL, XWORK2_FULL
48 !
49 USE MODD_PGD_GRID,       ONLY : NL, XGRID_PAR, NGRID_PAR, XMESHLENGTH
50 USE MODN_PGD_GRID
51 USE MODD_SURF_ATM_GRID_n, ONLY : XLAT, XLON, XMESH_SIZE, XJPDIR
52 USE MODD_SURF_ATM_n,     ONLY : NDIM_FULL, NSIZE_FULL
53 USE MODD_CSTS,           ONLY : XPI, XRADIUS
54 !
55 USE MODI_DEFAULT_GRID
56 USE MODI_GRID_FROM_FILE
57 USE MODI_OPEN_NAMELIST
58 USE MODI_TEST_NAM_VAR_SURF
59 USE MODI_CLOSE_NAMELIST
60 USE MODI_GET_LUOUT
61 USE MODI_READ_NAM_GRIDTYPE
62 USE MODI_LATLON_GRID
63 !
64 USE MODE_POS_SURF
65 !
66 !
67 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
68 USE PARKIND1  ,ONLY : JPRB
69 !
70 USE MODI_ABOR1_SFX
71 !
72 USE MODI_PGD_GRID_IO_INIT
73 #ifdef MNH_PARALLEL
74 USE MODE_TOOLS_ll, ONLY : GET_MEAN_OF_COORD_SQRT_ll
75 !
76 USE MODI_GET_SIZE_FULL_n
77 USE MODI_SPLIT_GRID
78 USE MODD_CONF, ONLY : CPROGRAM
79 #endif
80 !
81 IMPLICIT NONE
82 !
83 !*    0.1    Declaration of dummy arguments
84 !            ------------------------------
85 !
86  CHARACTER(LEN=6),  INTENT(IN)   :: HPROGRAM   ! program calling the surface
87  CHARACTER(LEN=28), INTENT(IN)   :: HFILE      ! atmospheric file name
88  CHARACTER(LEN=6),  INTENT(IN)   :: HFILETYPE  ! atmospheric file type
89 LOGICAL,           INTENT(IN)   :: OGRID      ! .true. if grid is imposed by atm. model
90  CHARACTER(LEN=10), INTENT(OUT)  :: HGRID      ! grid type
91 INTEGER,           INTENT(OUT)  :: KGRID_PAR  ! size of PGRID_PAR
92 REAL, DIMENSION(:), POINTER     :: PGRID_PAR  ! parameters defining this grid
93 !
94 !
95 !*    0.2    Declaration of local variables
96 !            ------------------------------
97 !
98 INTEGER           :: ILUOUT ! output listing logical unit
99 INTEGER           :: ILUNAM ! namelist file  logical unit
100 LOGICAL           :: GFOUND ! Flag true if namelist is present
101 REAL(KIND=JPRB) :: ZHOOK_HANDLE
102 INTEGER :: IIMAX_ll, IJMAX_ll ! global size of son model
103 LOGICAL :: GRECT
104 !
105  INTEGER :: IXOR = 1            ! position of modified bottom left point
106  INTEGER :: IYOR = 1            ! according to initial grid
107  INTEGER :: IXSIZE = -999       ! number of grid meshes in initial grid to be
108  INTEGER :: IYSIZE = -999       ! covered by the modified grid
109  INTEGER :: IDXRATIO = 1        ! resolution ratio between modified grid
110  INTEGER :: IDYRATIO = 1        ! and initial grid
111 NAMELIST/NAM_INIFILE_CONF_PROJ/IXOR,IYOR,IXSIZE,IYSIZE,IDXRATIO,IDYRATIO
112 !
113 !*    0.3    Declaration of namelists
114 !            ------------------------
115 !
116 !------------------------------------------------------------------------------
117 !
118 !*    1.      Defaults
119 !             --------
120 !
121 IF (LHOOK) CALL DR_HOOK('PGD_GRID',0,ZHOOK_HANDLE)
122  CALL DEFAULT_GRID(HPROGRAM,CGRID)
123 !
124 YINIFILE  = '                         '
125 YFILETYPE = '      '
126 !
127 IF (OGRID) THEN
128   YINIFILE  = HFILE
129   YFILETYPE = HFILETYPE
130 END IF
131 !
132  CALL GET_LUOUT(HPROGRAM,ILUOUT)
133 !------------------------------------------------------------------------------
134 !
135 !*    2.      Open namelist
136 !             -------------
137 !
138 IF (.NOT. OGRID) THEN
139   CALL OPEN_NAMELIST(HPROGRAM,ILUNAM)
140 !
141 !------------------------------------------------------------------------------
142 !
143 !*    3.      Read grid type
144 !             --------------
145 !
146   CALL POSNAM(ILUNAM,'NAM_PGD_GRID',GFOUND,ILUOUT)
147   IF (GFOUND) READ(UNIT=ILUNAM,NML=NAM_PGD_GRID)
148 !
149 !------------------------------------------------------------------------------
150 !
151 !*    5.      Close namelist
152 !             --------------
153 !
154   CALL CLOSE_NAMELIST(HPROGRAM,ILUNAM)
155 !
156 END IF
157 !-------------------------------------------------------------------------------
158 !
159 !*       4.    check of grid and input file types
160 !              ----------------------------------
161
162  CALL TEST_NAM_VAR_SURF(ILUOUT,'CGRID',CGRID,'CONF PROJ ','NONE      ','LONLAT REG','CARTESIAN ','GAUSS     ',&
163           'IGN       ','LONLATVAL ')  
164  CALL TEST_NAM_VAR_SURF(ILUOUT,'YFILETYPE',YFILETYPE,'      ','MESONH','LFI   ','ASCII ')
165 !
166 !
167 !------------------------------------------------------------------------------
168 !
169 !*    5.      Initializes grid characteristics
170 !             --------------------------------
171 !
172 !*    5.1     From another file
173 !             -----------------
174 !
175 IF (LEN_TRIM(YFILETYPE)>0 .AND. LEN_TRIM(YINIFILE)>0 ) THEN
176   IF (YFILETYPE=='MESONH' .OR. YFILETYPE=='LFI   ' .OR. YFILETYPE=='ASCII ') THEN
177     CALL GRID_FROM_FILE(HPROGRAM,YINIFILE,YFILETYPE,OGRID,CGRID,NGRID_PAR,XGRID_PAR,NL)
178     HGRID     = CGRID
179     IF ( HGRID == "IGN       " .OR. HGRID == "GAUSS     " .OR. HGRID == "NONE      " ) THEN
180       GRECT = .FALSE.
181     ELSE
182       GRECT = .TRUE.
183     ENDIF
184     ! on lit la taille globale du modele fils dans la namelist
185     CALL OPEN_NAMELIST(HPROGRAM,ILUNAM)
186     CALL POSNAM(ILUNAM,'NAM_INIFILE_CONF_PROJ',GFOUND,ILUOUT)
187     IF (GFOUND) THEN 
188       READ(UNIT=ILUNAM,NML=NAM_INIFILE_CONF_PROJ)
189       IIMAX_ll = IXSIZE*IDXRATIO
190       IJMAX_ll = IYSIZE*IDYRATIO
191     ENDIF
192     CALL CLOSE_NAMELIST(HPROGRAM,ILUNAM)
193
194     !*    3.      Additional actions for I/O
195     !
196     IF (GFOUND) THEN 
197 #ifdef MNH_PARALLEL
198       CALL PGD_GRID_IO_INIT(HPROGRAM,NGRID_PAR,XGRID_PAR, HGRID, GRECT, IIMAX_ll, IJMAX_ll, IDXRATIO, IDYRATIO)
199 #else
200       CALL PGD_GRID_IO_INIT(HPROGRAM)
201 #endif
202       NDIM_FULL = NL
203     ELSE
204 #ifdef MNH_PARALLEL
205       CALL PGD_GRID_IO_INIT(HPROGRAM,NGRID_PAR,XGRID_PAR, HGRID, GRECT)
206 #else
207       CALL PGD_GRID_IO_INIT(HPROGRAM)
208 #endif
209     ENDIF
210     NSIZE     = NDIM_FULL
211 #ifdef MNH_PARALLEL
212     CALL GET_SIZE_FULL_n(HPROGRAM,NDIM_FULL,NSIZE_FULL)
213     NL = NSIZE_FULL
214 #else
215     NSIZE_FULL = NL
216 #endif
217   ELSE
218     CALL ABOR1_SFX('PGD_GRID: FILE TYPE NOT SUPPORTED '//HFILETYPE//' FOR FILE '//HFILE)
219   END IF
220   !we don't need to call SPLIT_GRID, the grid has been splitted in GRID_FROM_FILE
221 !
222 ELSE
223 !
224 !*    5.2     Grid not initialized
225 !             --------------------
226 !
227   IF (CGRID=='NONE      ' .OR. CGRID=='          ') THEN
228     CALL ABOR1_SFX('PGD_GRID: GRID TYPE NOT INITIALIZED, CGRID='//CGRID)
229
230 !
231 !*    5.3     Grid initialized
232 !             ----------------
233 !
234   ELSE
235 !
236     CALL READ_NAM_GRIDTYPE(HPROGRAM,CGRID,NGRID_PAR,XGRID_PAR,NL)
237     HGRID     = CGRID
238     !*    3.      Additional actions for I/O
239     !
240 #ifdef MNH_PARALLEL
241     CALL PGD_GRID_IO_INIT(HPROGRAM,NGRID_PAR,XGRID_PAR)
242 #else
243     CALL PGD_GRID_IO_INIT(HPROGRAM)
244 #endif
245     NDIM_FULL = NL
246     NSIZE     = NDIM_FULL
247 #ifdef MNH_PARALLEL
248     CALL GET_SIZE_FULL_n(HPROGRAM,NDIM_FULL,NSIZE_FULL)
249     NL = NSIZE_FULL
250 #else
251     NSIZE_FULL = NL
252 #endif
253   END IF
254 #ifdef MNH_PARALLEL
255   ! IF we are in PREP_PGD, we need to split the grid. Otherwise, the grid was read in parallel and is already splitted
256   IF ( CPROGRAM == 'PGD   ') THEN
257     CALL SPLIT_GRID('MESONH',NGRID_PAR,XGRID_PAR)
258   ENDIF
259 #endif
260
261 END IF
262 !
263 !
264 IF (.NOT.ALLOCATED(NINDEX)) THEN
265   ALLOCATE(NINDEX(NDIM_FULL))
266   NINDEX(:) = 0
267 ENDIF
268 NINDX2    = NDIM_FULL
269 ALLOCATE(NWORK(NDIM_FULL))
270 ALLOCATE(XWORK(NDIM_FULL))
271 ALLOCATE(XWORK2(NDIM_FULL,10))
272 ALLOCATE(XWORK3(NDIM_FULL,10,10))
273 IF (NRANK==NPIO) THEN
274   ALLOCATE(NWORK_FULL(NDIM_FULL))
275   ALLOCATE(XWORK_FULL(NDIM_FULL))
276   ALLOCATE(XWORK2_FULL(NDIM_FULL,10))
277 ELSE
278   ALLOCATE(NWORK_FULL(0))
279   ALLOCATE(XWORK_FULL(0))
280   ALLOCATE(XWORK2_FULL(0,0))
281 ENDIF
282 !
283 KGRID_PAR = NGRID_PAR
284 ALLOCATE(PGRID_PAR(KGRID_PAR))
285 PGRID_PAR = XGRID_PAR
286 !
287 !------------------------------------------------------------------------------
288 !
289 !*    6.      Latitude and longitude
290 !             ----------------------
291 !
292 ALLOCATE(XLAT       (NSIZE_FULL))
293 ALLOCATE(XLON       (NSIZE_FULL))
294 ALLOCATE(XMESH_SIZE (NSIZE_FULL))
295 ALLOCATE(XJPDIR     (NSIZE_FULL))
296  CALL LATLON_GRID(CGRID,NGRID_PAR,NSIZE_FULL,ILUOUT,XGRID_PAR,XLAT,XLON,XMESH_SIZE,XJPDIR)
297 !
298 !------------------------------------------------------------------------------
299 !
300 !*    7.      Average grid length (in degrees)
301 !             --------------------------------
302 !
303 !* in meters
304 #ifdef MNH_PARALLEL
305 CALL GET_MEAN_OF_COORD_SQRT_ll(XMESH_SIZE,NSIZE_FULL,NDIM_FULL,XMESHLENGTH)
306 #else
307 XMESHLENGTH = SUM ( SQRT(XMESH_SIZE) ) / NL
308 #endif
309 !
310 !* in degrees (of latitude)
311 XMESHLENGTH = XMESHLENGTH *180. / XPI / XRADIUS
312 IF (LHOOK) CALL DR_HOOK('PGD_GRID',1,ZHOOK_HANDLE)
313 !
314 !-------------------------------------------------------------------------------
315 !
316 END SUBROUTINE PGD_GRID