source: TOOLS/MOZAIC/src/MOZAIC/lmdz.f90 @ 3326

Last change on this file since 3326 was 3326, checked in by omamce, 7 years ago

O.M. : Utility to generate interpolatio weights for OASIS-MCT
File size: 26.9 KB
Line 
1! -*- Mode: f90 -*-
2PROGRAM lmdgrille
3!!!
4!!! Ecris les grilles LMDZ au format OASIS
5   !!
6   !! Attention :
7   !!  pour OASIS 3, le nord est en haut (j=jpai) et le sud en bas (j=1)
8   !!  pour OASI3-MCT, il faut conserver le sens du modèle
9   !!
10   !! Unite IEEE : 32
11   !!
12   USE declare
13   USE dimensions
14   USE mod_prih
15   USE fliocom
16   USE getincom
17   USE errioipsl
18   USE mod_inipar
19   IMPLICIT none
20   !!
21   INTEGER (kind=il) :: jai, jaj, jk, jc
22   REAL (kind=rl), DIMENSION (:, :), ALLOCATABLE :: xlont, xlatt, xlonu, xlatu, xlonv, xlatv, xlonf, xlatf
23   REAL (kind=rl), DIMENSION (:, :), ALLOCATABLE :: xlat_o2a
24   REAL (kind=rl), DIMENSION (:, :), ALLOCATABLE :: xmas, xare, xtmp, xfra
25   REAL (kind=rl), DIMENSION (:, :, :), ALLOCATABLE :: xclon, xclat
26   REAL (kind=rl), DIMENSION (:, :), ALLOCATABLE :: xaw, xae, yan, yas
27   INTEGER (kind=il) :: id_restart, id_restartphy, id_hist, id_grid, id_coo, no = 32, narg = 0, iargc
28   REAL (kind=rl) :: ra, zzz, xdel, ydel
29   CHARACTER (len=180) :: cfname, clmd, c_restart, c_restartphy, c_hist
30   CHARACTER (len=80) :: cmodel
31   CHARACTER (len=8) :: clon, clat, cmsk, csrf, clarg
32   REAL (kind=rl), DIMENSION (:, :), ALLOCATABLE :: flux_iceberg
33   INTEGER :: ia_sk=1_il, ja_sk=1_il
34   INTEGER :: nlmd
35   REAL (kind=rl), DIMENSION (:), ALLOCATABLE :: tab1d
36   !!
37   INTEGER (kind=il) :: il_grid, il_mask, il_area, il_frac, id
38   INTEGER (kind=il) :: ierr
39   !
40   LOGICAL :: l_build   = .FALSE.  !< Construction a partir de rien
41   LOGICAL :: l_restart = .FALSE.  !< Construction a partir de restart et restartphy
42   LOGICAL :: l_read    = .FALSE.  !< Construction a partir de grids.nc, areas.nc, masks.nc
43   !LOGICAL :: lec_msk = .TRUE., lec_srf = .TRUE., lec_coo = .TRUE., lec_fra = .TRUE.
44   LOGICAL :: la_retourn_lat !< Retournement nor/sud des champs lmdz
45   CHARACTER (LEN=1) :: cl_rk
46   !!
47   CALL inipar
48   
49   !!
50   IF ( jpai .GT. 100_il ) ia_sk =  2
51   IF ( jpai .GT. 200_il ) ia_sk =  5
52   IF ( jpai .GT. 500_il ) ia_sk = 10
53
54   IF ( jpaj .GT. 100_il ) ia_sk =  2
55   IF ( jpaj .GT. 200_il ) ja_sk =  5
56   IF ( jpaj .GT. 500_il ) ja_sk = 10
57   !!
58   ALLOCATE (xlont (jpai , jpaj ), STAT=ierr) ; CALL chk_allo (ierr, 'xlont', lreset=.TRUE., crout='lmdz')
59   ALLOCATE (xlatt (jpai , jpaj ), STAT=ierr) ; CALL chk_allo (ierr, 'xlatt')
60   ALLOCATE (xlonu (jpaiu, jpaju), STAT=ierr) ; CALL chk_allo (ierr, 'xlonu')
61   ALLOCATE (xlatu (jpaiu, jpaju), STAT=ierr) ; CALL chk_allo (ierr, 'xlatu')
62   ALLOCATE (xlonv (jpaiv, jpajv), STAT=ierr) ; CALL chk_allo (ierr, 'xlonv')
63   ALLOCATE (xlatv (jpaiv, jpajv), STAT=ierr) ; CALL chk_allo (ierr, 'xlatv')
64   ALLOCATE (xmas  (jpai , jpaj ), STAT=ierr)  ; CALL chk_allo (ierr, 'xmas')
65   ALLOCATE (xare  (jpai , jpaj ), STAT=ierr)  ; CALL chk_allo (ierr, 'xare')
66   ALLOCATE (xfra  (jpai , jpaj ), STAT=ierr)  ; CALL chk_allo (ierr, 'xfra')
67   ALLOCATE (xtmp  (jpai , jpaj )) ; CALL chk_allo (ierr, 'xtmp')
68   ALLOCATE (xclon (jpai , jpaj, 4), STAT=ierr) ; CALL chk_allo (ierr, 'xclon')
69   ALLOCATE (xclat (jpai , jpaj, 4), STAT=ierr) ; CALL chk_allo (ierr, 'xclat')
70   ALLOCATE (xaw   (jpai , jpaj), STAT=ierr) ; CALL chk_allo (ierr, 'xaw')
71   ALLOCATE (xae   (jpai , jpaj), STAT=ierr) ; CALL chk_allo (ierr, 'xae')
72   ALLOCATE (yan   (jpai , jpaj), STAT=ierr) ; CALL chk_allo (ierr, 'yan')
73   ALLOCATE (yas   (jpai , jpaj), STAT=ierr) ; CALL chk_allo (ierr, 'yas')
74   ALLOCATE (flux_iceberg (jpai, jpaj)) ; CALL chk_allo (ierr, 'flux_iceberg')
75
76   !!
77   ra =  r_earth
78   !!
79   clon = 'nav_lon' ; clat = 'nav_lat' ; cmsk = 'msks' ; csrf = 'aire' ; clmd =  'grillelmd.nc'
80   !!
81   !CALL ipsldbg (.TRUE.)
82   !!
83   WRITE (UNIT=cmodel, FMT='("LMDZ ", 1I3.3, "x", 1I3.3 )') jpai, jpaj
84   WRITE ( nout, *) cmodel
85   !!
86   !! Lecture de la ligne de commande
87   narg = iargc()
88   IF ( narg > 0 ) THEN
89      CALL getarg ( 1, clarg)
90      SELECT CASE (TRIM(clarg))
91      CASE ('-b')
92         l_build = .TRUE.
93   
94      CASE ('-l')
95         l_restart = .TRUE. 
96         c_restart    = 'restart.nc'
97         c_restartphy = 'restartphy.nc'
98!-$$         c_hist       = 'histmth.nc'
99!-$$         csrf         = 'AIRE'
100         c_hist       = 'gridarea_zoom_correct.nc'
101         csrf         = 'area'
102         WRITE ( nout, *) 'Lecture LMDZ dans restart : ', TRIM (c_restart)
103
104      CASE ('-r') 
105         l_read  = .TRUE. 
106         clmd = 'o2a.diag.nc' 
107         cmsk = 'OceMask' ; csrf = 'Surface' 
108         WRITE ( nout, *) 'Relecture du masque calcule par mozaic'
109
110      CASE Default
111         IF ( INDEX ( TRIM(clarg), '.def' ) /= 0 ) THEN
112            WRITE (nout,*) 'Lecture des parametres dans ', TRIM(clarg)
113            CALL getin_name ( TRIM (clarg) )
114         ELSE
115            WRITE (nout,*) 'Lecture des parametres dans run.def'
116         ENDIF
117      END SELECT
118   END IF
119   !!
120   WRITE ( nout,*) 'Dimension atmopshere : ', jpai, jpaj, jpan
121   WRITE ( nout,*) 'Dimension ocean      : ', jpoi, jpoj, jpon
122!-$$   WRITE ( nout,*) 'Lecture mask  : ', lec_msk
123!-$$   WRITE ( nout,*) 'Lecture coord : ', lec_coo
124   !!
125   !!
126   WRITE ( unit = cl_rk, fmt = '(1I1)' ) rk_in
127   cfname = 'lmd.grids' // TRIM (c_suffix) // '.r' // cl_rk
128   OPEN ( unit = no, file = TRIM(cfname), form = 'unformatted', action = 'write', position = 'rewind', status = 'unknown' )
129   !
130   !CALL ipsldbg ( .TRUE. )
131   !!
132   IF ( l_build ) THEN
133      WRITE (nout,*) 'Construction de la grille ex nihilo'
134      ! Lon, lat - regulier uniquement
135      ! Contruction avec le Sud en bas de la carte
136      DO jai = 1, jpai
137         xlont (jai,:) = -180.0_rl + (REAL (jai-1, rl)        ) / REAL (jpai  , rl) * 360.0_rl
138         xaw   (jai,:) = -180.0_rl + (REAL (jai-1, rl)-0.5_rl ) / REAL (jpai  , rl) * 360.0_rl
139         xae   (jai,:) = -180.0_rl + (REAL (jai-1, rl)+0.5_rl ) / REAL (jpai  , rl) * 360.0_rl
140      END DO
141     
142      DO jaj = 1, jpaj
143         xlatt (:,jaj) =   90.0_rl - ( REAL (jaj-1, rl)        ) / REAL (jpaj-1, rl) * 180.0_rl
144         yas   (:,jaj) =   90.0_rl - ( REAL (jaj-1, rl)+0.5_rl ) / REAL (jpaj-1, rl) * 180.0_rl
145         yan   (:,jaj) =   90.0_rl - ( REAL (jaj-1, rl)-0.5_rl ) / REAL (jpaj-1, rl) * 180.0_rl
146      END DO
147      !
148      yas = MIN ( 90.0_rl, MAX (-90.0_rl, yas))
149      yan = MIN ( 90.0_rl, MAX (-90.0_rl, yan))
150      !
151      xclon (:,:,1) = xae
152      xclon (:,:,2) = xaw
153      xclon (:,:,3) = xaw
154      xclon (:,:,4) = xae
155      !
156      xclat (:,:,1) = yas
157      xclat (:,:,2) = yas
158      xclat (:,:,3) = yan
159      xclat (:,:,4) = yan
160     
161      DO jk = 1, 4
162         xclon(:,:,jk) = clo_lon ( xclon(:,:,jk), xlont(:,:) )
163      END DO
164     
165      ! Masque et fractions
166      xmas = 0.0_rl
167      xfra = 1.0_rl
168
169      ! Surfaces
170      CALL surface ( 'globe' )
171
172   END IF
173   
174   IF ( l_restart ) THEN
175      !! Lecture restart LMDZ : le sud est en bas
176      nlmd =  jpait*(jpajv-1) + 2
177      WRITE (nout,*) 'Lecture restart LMD '
178      WRITE (nout,*) 'Dimension 1D : ', nlmd
179      ALLOCATE ( tab1d (nlmd)) ; CALL chk_allo (ierr, 'tab1d')
180      CALL flioopfd ( TRIM(c_restart)   , id_restart)
181      CALL flioopfd ( TRIM(c_restartphy), id_restartphy)
182      !
183      WRITE (nout,*) 'Lecture longitude ', jpait*(jpajv-1), (/jpait, jpajv-1/)
184      CALL fliogetv (id_restartphy, 'longitude', tab1d (:))
185      xlont (:,2:jpai-1) = RESHAPE ( tab1d (2:nlmd-1), (/jpait, jpajv-1/) )
186      xlont (:, 1  ) = xlont (:,jpaj/2)
187      xlont (:,jpaj) = xlont (:,jpaj/2)
188      WRITE (nout,*) 'Lecture latitude ', jpait*(jpajv-1), (/jpait, jpajv-1/) 
189      CALL fliogetv (id_restartphy, 'latitude' , tab1d (:))
190      xlatt (:,2:jpai-1) = RESHAPE ( tab1d (2:nlmd-1), (/jpait, jpajv-1/) )
191      xlatt (:,1) = 90.0_rl ; xlatt (:,jpaj) = -90.0_rl
192
193      WRITE (nout,*) 'Lecture FOCE ', jpait*(jpajv-1), (/jpait, jpajv-1/)
194      CALL fliogetv (id_restartphy, 'FOCE' , tab1d (:))
195      xfra (:,2:jpai-1) = RESHAPE ( tab1d (1:nlmd-1), (/jpait, jpajv-1/) )
196      WRITE (nout,*) 'Lecture FSIC ', jpait*(jpajv-1), (/jpait, jpajv-1/)
197      CALL fliogetv (id_restartphy, 'FSIC' , tab1d (:))
198      xfra (:,2:jpai-1) = xfra (:,2:jpai-1) & 
199         &              + RESHAPE ( tab1d (1:nlmd-1), (/jpait, jpajv-1/) )
200      xfra (:, 1   ) = 1.0_rl
201      xfra (:, jpaj) = 0.0_rl
202
203      xmas (:,:) = 1.0_rl
204      WHERE ( xfra (:,:) .GT. 0.0_rl) xmas (:,:) = 0.0_rl
205
206      CALL fliogetv (id_restart, 'rlonu', xlonu(:,1) )
207      CALL fliogetv (id_restart, 'rlatu', xlatu(1,:) )
208      CALL fliogetv (id_restart, 'rlonv', xlonv(:,1) )
209      CALL fliogetv (id_restart, 'rlatv', xlatv(1,:) )
210
211      ! Eventuellement conversion
212      IF ( MAXVAL(ABS(xlonu)) .LT. 10.0_rl ) THEN
213         xlonu (:,:) = SPREAD ( xlonu(:,1), DIM=2, ncopies=jpaju) * dar
214         xlatu (:,:) = SPREAD ( xlatu(1,:), DIM=1, ncopies=jpaiu) * dar
215         xlonv (:,:) = SPREAD ( xlonv(:,1), DIM=2, ncopies=jpajv) * dar
216         xlatv (:,:) = SPREAD ( xlatv(1,:), DIM=1, ncopies=jpaiv) * dar
217      END IF
218
219      !
220      xaw  (2:jpai  ,:) = xlonu (1:jpai-1,1:jpaj)
221      xaw  (1       ,:) = xlonu (jpai    ,1:jpaj) - 360.0_rl
222
223      xae  (1:jpai  ,:) = xlonu (1:jpai  ,1:jpaj)
224
225      !
226      yan (1:jpai,2:jpaj  ) = xlatv(1:jpai,1:jpaj-1)
227      yan ( :    ,  1     ) = 90.0_rl
228
229      yas (1:jpai,1:jpaj-1) = xlatv (1:jpai, 1:jpaj-1)
230      yas ( :    ,jpaj    ) = -90.0_rl
231
232      xclon (:,:,1) = xae
233      xclon (:,:,2) = xaw
234      xclon (:,:,3) = xaw
235      xclon (:,:,4) = xae
236      !
237      xclat (:,:,1) = yas
238      xclat (:,:,2) = yas
239      xclat (:,:,3) = yan
240      xclat (:,:,4) = yan
241
242      ! Surfaces
243      WRITE (nout,*) 'Lecture hist LMD ', TRIM(c_hist)
244      CALL flioopfd ( TRIM(c_hist), id_hist)
245      !
246      CALL fliogetv (id_hist, TRIM(csrf), xare(:,:), start=(/1,1,1/), count=(/jpai,jpaj,1/) )
247     
248      !! Bug surface dans le zoom de Yves et Rng
249      !CALL surface ( 'pole' )
250     
251   ENDIF
252
253   IF ( l_read ) THEN
254      ! Le sud n'est pas forcement au bon endroit
255      CALL flioopfd ('grids' // TRIM(c_suffix) // '.nc', id)
256      CALL fliogetv (id, 'tlmd.lon', xlont)
257      CALL fliogetv (id, 'tlmd.lat', xlatt)
258      CALL fliogetv (id, 'tlmd.clo', xclon)
259      CALL fliogetv (id, 'tlmd.cla', xclat)
260      CALL flioclo (id)
261     
262      CALL flioopfd ('areas' // TRIM(c_suffix) // '.nc', id)
263      CALL fliogetv (id, 'tlmd.srf', xare)
264      CALL flioclo (id)
265     
266      WRITE ( nout, *) 'Lecture des fractions dans ' // TRIM(clmd) // ' variable ' // TRIM(cmsk) 
267      ! Les fractions ne sont pas forcément dans le même sens que les autres fichiers ... !
268      CALL flioopfd ( TRIM (clmd), id)
269      CALL fliogetv ( id, TRIM (cmsk), xfra)
270      WRITE ( UNIT = nout, fmt = *) 'Msk lmd ', TRIM ( cmsk)
271      CALL prihre ( xfra, 1.0_rl, nout)
272      IF ( TRIM (cmsk) == 'phis' ) THEN
273         xtmp = xmas
274         xmas = 0.0_rl
275         WHERE ( xtmp .GT. 0.1_rl ) xmas = 1.0_rl
276         xfra = xtmp
277      END IF
278      IF ( TRIM (cmsk) == 'OceMask' ) THEN
279         xmas = 0.0_rl
280         WHERE ( xfra .EQ. 0.0_rl ) xmas = 1.0_rl
281         ALLOCATE ( xlat_o2a(jpai,jpaj))
282         CALL fliogetv ( id, 'nav_lat', xlat_o2a)
283         !
284         IF (     xlat_o2a(jpai/2,1) .LT.  xlat_o2a(jpai/2,jpaj) .AND. xlatt(jpai/2,1) .GT. xlatt(jpai/2,jpaj) &
285            .OR.  xlat_o2a(jpai/2,1) .GT.  xlat_o2a(jpai/2,jpaj) .AND. xlatt(jpai/2,1) .LT. xlatt(jpai/2,jpaj) ) THEN
286            xfra ( :,:) = xfra (:, jpaj:1:-1) 
287            xmas ( :,:) = xmas (:, jpaj:1:-1) 
288         END IF
289      END IF
290   END IF
291   !!
292   !! Inversion eventuelle des indices J
293   !!
294   IF ( la_nortop ) THEN
295      WRITE (unit = nout, fmt = *) "Cas nord en haut"
296      IF ( xlatt(jpai/2,1) .GT. xlatt(jpai/2,jpaj) ) THEN
297         WRITE (nout,*) "Il faut retourner les indices J"
298         la_retourn_lat  = .TRUE.
299      ELSE
300         WRITE (nout,*) "Indices J dans le bon ordre"
301         la_retourn_lat  = .FALSE.
302      END IF
303   ELSE
304      WRITE (unit = nout, fmt = *) "Cas nord en bas"
305      IF ( xlatt(jpai/2,1) .LT. xlatt(jpai/2,jpaj) ) THEN
306         WRITE (nout,*) "Il faut retourner les indices J"
307         la_retourn_lat  = .TRUE.
308      ELSE
309         WRITE (nout,*) "Indices J dans le bon ordre"
310         la_retourn_lat  = .FALSE.
311      END IF
312   END IF
313
314   IF ( la_retourn_lat ) THEN
315      WRITE (nout,*) "Retournement des indices J"
316      xlont (:,:)   =  xlont (1:jpai , jpaj :1:-1)
317      xlonu (:,:)   =  xlonu (1:jpaiu, jpaju:1:-1)
318      xlonv (:,:)   =  xlonv (1:jpaiv, jpajv:1:-1)
319      xlatt (:,:)   =  xlatt (1:jpai , jpaj :1:-1)
320      xlatu (:,:)   =  xlatu (1:jpaiu, jpaju:1:-1)
321      xlatv (:,:)   =  xlatv (1:jpaiv, jpajv:1:-1)
322      xmas  (:,:)   =  xmas  (1:jpai , jpaj :1:-1)
323      xfra  (:,:)   =  xfra  (1:jpai , jpaj :1:-1)
324      xare  (:,:)   =  xare  (1:jpai , jpaj :1:-1)
325      xclon (:,:,:) =  xclon (1:jpai , jpaj :1:-1,:)
326      xclat (:,:,:) =  xclat (1:jpai , jpaj :1:-1,:)
327
328      xaw   (:,:)   =  xaw   (1:jpai , jpaj :1:-1)
329      xae   (:,:)   =  xae   (1:jpai , jpaj :1:-1)
330      yas   (:,:)   =  yas   (1:jpai , jpaj :1:-1)
331      yan   (:,:)   =  yan   (1:jpai , jpaj :1:-1)
332
333      xclon (:,:,1) = xaw (:,:) 
334      xclon (:,:,2) = xae (:,:) 
335      xclon (:,:,3) = xae (:,:) 
336      xclon (:,:,4) = xaw (:,:) 
337      !
338      xclat (:,:,1) = yan (:,:) 
339      xclat (:,:,2) = yan (:,:) 
340      xclat (:,:,3) = yas (:,:) 
341      xclat (:,:,4) = yas (:,:) 
342
343   END IF
344
345   WRITE ( unit = nout, fmt = *) 'Longitudes T'
346   CALL prihre ( xlont, 1._rl, nout)
347
348   WRITE ( unit = no) 't'//TRIM(camod)//'.lon'
349   WRITE ( unit = no) REAL ( xlont(1:jpai,1:jpaj), kind = rk_in )
350   !
351   WRITE ( unit = nout, fmt = *) 'Latitudes T'
352   CALL prihre ( xlatt, 1._rl, nout)
353
354   WRITE ( unit = no) 't'//TRIM(camod)//'.lat'
355   WRITE ( unit = no) REAL (xlatt(1:jpai,1:jpaj), kind = rk_in )
356   !
357   WRITE ( unit = nout, fmt = *) 'Corners'
358   WRITE ( unit = no) 't'//TRIM(camod)//'.clo'
359   WRITE ( unit = no) REAL (xclon, kind = rk_in )
360   WRITE ( unit = no) 't'//TRIM(camod)//'.cla'
361   WRITE ( unit = no) REAL (xclat, kind = rk_in )
362
363   DO jc = 1, 4
364      WRITE ( unit = nout, fmt = *) 'Coins - Longitudes ', jc
365      CALL prihre ( xclon(:,:,jc), 1._rl, nout)
366      WRITE ( unit = nout, fmt = *) 'Coins - Latitudes ', jc
367      CALL prihre ( xclat(:,:,jc), 1._rl, nout)
368   END DO
369
370   WRITE ( unit = nout, fmt = *) 'Fractions (%)'
371   CALL prihre ( xfra, 100.0_rl, nout)
372
373   WRITE ( unit = nout, fmt = *) 'Surface '
374   CALL prihre (xare, 0.0_rl, nout)
375   zzz = SUM ( xare) / ( ra * ra * rpi * 4.0_rl )
376   WRITE (UNIT = nout, FMT = *) 'Surface normalisee : ', zzz
377
378   !!
379   !! Ecriture OASIS NetCDF
380   CALL fliocrfd ( 'lmd.grids' // TRIM(c_suffix), (/'jpia   ', 'jpja   ', 'corners'/), (/jpai, jpaj, 4/), il_grid, mode='REPLACE')
381   CALL flioputa ( il_grid, '?', 'Model', cmodel)
382   CALL flioputa ( il_grid, '?', 'Comment', TRIM(c_comment) )
383
384   CALL fliodefv ( il_grid, 'tlmd.lon', (/1, 2/), v_t=flio_r8 )
385   CALL flioputa ( il_grid, 'tlmd.lon', 'units'     , 'degrees_east' )
386   CALL flioputa ( il_grid, 'tlmd.lon', 'title'     , 'lmdz-t longitudes' )
387   CALL flioputa ( il_grid, 'tlmd.lon', 'grid_type' , 'P' )
388   CALL flioputa ( il_grid, 'tlmd.lon', 'overlap'   , 0_i_4 )
389
390   CALL fliodefv ( il_grid, 'tlmd.lat', (/1, 2/), v_t=flio_r8 )
391   CALL flioputa ( il_grid, 'tlmd.lat', 'units'    , 'degrees_north' )
392   CALL flioputa ( il_grid, 'tlmd.lat', 'title'     ,'lmdz-t latitudes' )
393   CALL flioputa ( il_grid, 'tlmd.lat', 'grid_type' ,'P' )
394   CALL flioputa ( il_grid, 'tlmd.lat', 'overlap'   , 0_i_4 )
395
396   CALL fliodefv ( il_grid, 'tlmd.clo', (/1, 2, 3/), v_t=flio_r8 )
397   CALL flioputa ( il_grid, 'tlmd.clo', 'units'    , 'degrees_east' )
398   CALL flioputa ( il_grid, 'tlmd.clo', 'title'     ,'Longitudes for t-cell corner' )
399   CALL flioputa ( il_grid, 'tlmd.clo', 'grid_type' ,'P' )
400   CALL flioputa ( il_grid, 'tlmd.clo', 'overlap'   , 0_i_4 )
401
402   CALL fliodefv ( il_grid, 'tlmd.cla', (/1, 2, 3/), v_t=flio_r8 )
403   CALL flioputa ( il_grid, 'tlmd.cla', 'units'    , 'degrees_north' )
404   CALL flioputa ( il_grid, 'tlmd.cla', 'title'     ,'Latitudes for t-cell corner' )
405   CALL flioputa ( il_grid, 'tlmd.cla', 'grid_type' ,'P' )
406   CALL flioputa ( il_grid, 'tlmd.cla', 'overlap'   , 0_i_4 )
407
408   CALL fliodefv ( il_grid, 'aone.lon', (/1, 2/), v_t=flio_r8 )
409   CALL flioputa ( il_grid, 'aone.lon', 'units'     , 'degrees_east' )
410   CALL flioputa ( il_grid, 'aone.lon', 'title'     , 'lmdz-t longitudes' )
411   CALL flioputa ( il_grid, 'aone.lon', 'grid_type' , 'P' )
412   CALL flioputa ( il_grid, 'aone.lon', 'overlap'   , 0_i_4 )
413
414   CALL fliodefv ( il_grid, 'aone.lat', (/1, 2/), v_t=flio_r8 )
415   CALL flioputa ( il_grid, 'aone.lat', 'units'    , 'degrees_north' )
416   CALL flioputa ( il_grid, 'aone.lat', 'title'     ,'lmdz-t latitudes' )
417   CALL flioputa ( il_grid, 'aone.lat', 'grid_type' ,'P' )
418   CALL flioputa ( il_grid, 'aone.lat', 'overlap'   , 0_i_4 )
419
420   CALL fliodefv ( il_grid, 'aone.clo', (/1, 2, 3/), v_t=flio_r8 )
421   CALL flioputa ( il_grid, 'aone.clo', 'units'    , 'degrees_east' )
422   CALL flioputa ( il_grid, 'aone.clo', 'title'     ,'Longitudes for t-cell corner' )
423   CALL flioputa ( il_grid, 'aone.clo', 'grid_type' ,'P' )
424   CALL flioputa ( il_grid, 'aone.clo', 'overlap'   , 0_i_4 )
425
426   CALL fliodefv ( il_grid, 'aone.cla', (/1, 2, 3/), v_t=flio_r8 )
427   CALL flioputa ( il_grid, 'aone.cla', 'units'    , 'degrees_north' )
428   CALL flioputa ( il_grid, 'aone.cla', 'title'     ,'Latitudes for t-cell corner' )
429   CALL flioputa ( il_grid, 'aone.cla', 'grid_type' ,'P' )
430   CALL flioputa ( il_grid, 'aone.cla', 'overlap'   , 0_i_4 )
431
432   CALL fliodefv ( il_grid, 'afra.lon', (/1, 2/), v_t=flio_r8 )
433   CALL flioputa ( il_grid, 'afra.lon', 'units'     , 'degrees_east' )
434   CALL flioputa ( il_grid, 'afra.lon', 'title'     , 'lmdz-t longitudes' )
435   CALL flioputa ( il_grid, 'afra.lon', 'grid_type' , 'P' )
436   CALL flioputa ( il_grid, 'afra.lon', 'overlap'   , 0_i_4 )
437
438   CALL fliodefv ( il_grid, 'afra.lat', (/1, 2/), v_t=flio_r8 )
439   CALL flioputa ( il_grid, 'afra.lat', 'units'    , 'degrees_north' )
440   CALL flioputa ( il_grid, 'afra.lat', 'title'     ,'lmdz-t latitudes' )
441   CALL flioputa ( il_grid, 'afra.lat', 'grid_type' ,'P' )
442   CALL flioputa ( il_grid, 'afra.lat', 'overlap'   , 0_i_4 )
443
444   CALL fliodefv ( il_grid, 'afra.clo', (/1, 2, 3/), v_t=flio_r8 )
445   CALL flioputa ( il_grid, 'afra.clo', 'units'    , 'degrees_east' )
446   CALL flioputa ( il_grid, 'afra.clo', 'title'     ,'Longitudes for t-cell corner' )
447   CALL flioputa ( il_grid, 'afra.clo', 'grid_type' ,'P' )
448   CALL flioputa ( il_grid, 'afra.clo', 'overlap'   , 0_i_4 )
449
450   CALL fliodefv ( il_grid, 'afra.cla', (/1, 2, 3/), v_t=flio_r8 )
451   CALL flioputa ( il_grid, 'afra.cla', 'units'    , 'degrees_north' )
452   CALL flioputa ( il_grid, 'afra.cla', 'title'     ,'Latitudes for t-cell corner' )
453   CALL flioputa ( il_grid, 'afra.cla', 'grid_type' ,'P' )
454   CALL flioputa ( il_grid, 'afra.cla', 'overlap'   , 0_i_4 )
455
456   CALL flioputv ( il_grid, 'tlmd.lon', xlont )
457   CALL flioputv ( il_grid, 'tlmd.lat', xlatt )
458   CALL flioputv ( il_grid, 'tlmd.clo', xclon )
459   CALL flioputv ( il_grid, 'tlmd.cla', xclat )
460   CALL flioputv ( il_grid, 'aone.lon', xlont )
461   CALL flioputv ( il_grid, 'aone.lat', xlatt )
462   CALL flioputv ( il_grid, 'aone.clo', xclon )
463   CALL flioputv ( il_grid, 'aone.cla', xclat )
464   CALL flioputv ( il_grid, 'afra.lon', xlont )
465   CALL flioputv ( il_grid, 'afra.lat', xlatt )
466   CALL flioputv ( il_grid, 'afra.clo', xclon )
467   CALL flioputv ( il_grid, 'afra.cla', xclat )
468   !
469   CALL flioclo ( il_grid)
470   !
471   CLOSE ( unit = no)
472   !!
473   WRITE ( unit = nout, fmt = *) 'Masque atm ', i_4
474   WRITE ( unit = cl_rk, fmt = '(1I1)' ) i_4
475   cfname = 'lmd.masks' // TRIM (c_suffix) // '.i' // cl_rk
476   OPEN ( unit = no, file = cfname, form = 'unformatted', action = 'write', position = 'rewind', status = 'unknown' )
477   CALL prihin ( INT (xmas,KIND=il), 1_il, nout)
478   WRITE ( unit = no) 't'//TRIM(camod)//'.msk'
479   WRITE ( unit = no) INT ( xmas, KIND = i_4)
480   WRITE ( unit = no) 'u'//TRIM(camod)//'.msk'
481   WRITE ( unit = no) INT ( xmas, KIND = i_4)
482   WRITE ( unit = no) 'v'//TRIM(camod)//'.msk'
483   WRITE ( unit = no) INT ( xmas, KIND = i_4)
484   CLOSE ( unit = no)
485   !!
486   WRITE ( unit = nout, fmt = *) 'Masque atm ', i_8
487   WRITE ( unit = cl_rk, fmt = '(1I1)' ) i_8
488   cfname = 'lmd.masks' // TRIM (c_suffix) // '.i' // cl_rk
489   OPEN ( unit = no, file = cfname, form = 'unformatted', action = 'write', position = 'rewind', status = 'unknown' )
490   CALL prihin ( INT (xmas,KIND=il), 1_il, nout)
491   WRITE ( unit = no) 't'//TRIM(camod)//'.msk'
492   WRITE ( unit = no) INT ( xmas, KIND = i_8)
493   WRITE ( unit = no) 'u'//TRIM(camod)//'.msk'
494   WRITE ( unit = no) INT ( xmas, KIND = i_8)
495   WRITE ( unit = no) 'v'//TRIM(camod)//'.msk'
496   WRITE ( unit = no) INT ( xmas, KIND = i_8)
497   CLOSE ( unit = no)
498   !!
499   CALL fliocrfd ( 'lmd.masks' // TRIM(c_suffix), (/'jpia', 'jpja' /), (/jpai, jpaj/), il_mask, mode='REPLACE' )
500   CALL flioputa ( il_mask, '?', 'Model', cmodel )
501   CALL flioputa ( il_mask, '?', 'Comment', TRIM(c_comment) )
502
503   CALL fliodefv ( il_mask, 'tlmd.msk', (/1, 2/), v_t=flio_i4 )
504   CALL flioputa ( il_mask, 'tlmd.msk', 'units'     , 'n' )
505   CALL flioputa ( il_mask, 'tlmd.msk', 'land_value', 1_i_4)
506   CALL flioputa ( il_mask, 'tlmd.msk', 'sea_value' , 0_i_4)
507   CALL flioputa ( il_mask, 'tlmd.msk', 'title'     , 'lmdz-t land-sea mask' )
508
509   CALL fliodefv ( il_mask, 'aone.msk', (/1, 2/), v_t=flio_i4 )
510   CALL flioputa ( il_mask, 'aone.msk', 'units'     , 'n' )
511   CALL flioputa ( il_mask, 'aone.msk', 'land_value', 1_i_4)
512   CALL flioputa ( il_mask, 'aone.msk', 'sea_value' , 0_i_4)
513   CALL flioputa ( il_mask, 'aone.msk', 'title'     , 'lmdz-t land-sea mask, forced to not masked' )
514
515   CALL fliodefv ( il_mask, 'afra.msk', (/1, 2/), v_t=flio_i4 )
516   CALL flioputa ( il_mask, 'afra.msk', 'units'     , 'n' )
517   CALL flioputa ( il_mask, 'afra.msk', 'land_value', 1_i_4)
518   CALL flioputa ( il_mask, 'afra.msk', 'sea_value' , 0_i_4)
519   CALL flioputa ( il_mask, 'afra.msk', 'title'     , 'lmdz-t land-sea mask' )
520
521   !!
522   CALL flioputv ( il_mask, 'tlmd.msk', xmas )
523   CALL flioputv ( il_mask, 'aone.msk', 0.0*xmas )
524   CALL flioputv ( il_mask, 'afra.msk', xmas )
525
526   CALL flioclo ( il_mask )
527   !!
528   WRITE ( unit = cl_rk, fmt = '(1I1)' ) rk_in
529   cfname = 'lmd.areas' // TRIM (c_suffix) // '.r' // cl_rk
530   OPEN ( unit = no, file = cfname, form = 'unformatted', action = 'write', position = 'rewind', status = 'unknown' )
531   !
532
533   WRITE (unit = no) 't'//TRIM(camod)//'.srf'
534   WRITE (unit = no) REAL ( xare, KIND = rk_in)
535   !
536   CLOSE ( unit = no)
537   !!
538   CALL fliocrfd ( 'lmd.areas' // TRIM(c_suffix), (/'jpia', 'jpja' /), (/jpai, jpaj/), il_area, mode='REPLACE' )
539   CALL flioputa ( il_area, '?', 'Model', cmodel )
540   CALL flioputa ( il_area, '?', 'Comment', TRIM(c_comment) )
541
542   CALL fliodefv ( il_area, 'tlmd.srf', (/1, 2/), v_t=flio_r8 )
543   CALL flioputa ( il_area, 'tlmd.srf', 'units'  , 'm^2' )
544   CALL flioputa ( il_area, 'tlmd.srf', 'title', 'lmdz-t mesh surface')
545
546   CALL fliodefv ( il_area, 'aone.srf', (/1, 2/), v_t=flio_r8 )
547   CALL flioputa ( il_area, 'aone.srf', 'units'  , 'm^2' )
548   CALL flioputa ( il_area, 'aone.srf', 'title', 'lmdz-t mesh surface, set to 1')
549
550   CALL fliodefv ( il_area, 'afra.srf', (/1, 2/), v_t=flio_r8 )
551   CALL flioputa ( il_area, 'afra.srf', 'units'  , 'm^2' )
552   CALL flioputa ( il_area, 'afra.srf', 'title', 'lmdz-t mesh surface, set to 1')
553
554   CALL flioputv ( il_area, 'tlmd.srf', xare )
555   CALL flioputv ( il_area, 'aone.srf', 1.0_rl + 0.0_rl*xare )
556   CALL flioputv ( il_area, 'afra.srf', xare*xfra )
557
558   CALL flioclo ( il_area)
559   !
560   !!
561   CALL fliocrfd ( 'lmd.fracs' // TRIM(c_suffix), (/'jpia', 'jpja' /), (/jpai, jpaj/), il_frac, mode='REPLACE' )
562   CALL flioputa ( il_frac, '?', 'Model', cmodel )
563   CALL flioputa ( il_frac, '?', 'Comment', TRIM(c_comment) )
564
565   CALL fliodefv ( il_frac, 'tlmd.fra', (/1, 2/), v_t=flio_r8 )
566   CALL flioputa ( il_frac, 'tlmd.fra', 'units'  , 'm^2' )
567   CALL flioputa ( il_frac, 'tlmd.fra', 'title', 'lmdz-t mesh fraction')
568
569   CALL fliodefv ( il_frac, 'aone.fra', (/1, 2/), v_t=flio_r8 )
570   CALL flioputa ( il_frac, 'aone.fra', 'units'  , 'm^2' )
571   CALL flioputa ( il_frac, 'aone.fra', 'title', 'lmdz-t mesh fraction, set to 1')
572
573   CALL fliodefv ( il_frac, 'afra.fra', (/1, 2/), v_t=flio_r8 )
574   CALL flioputa ( il_frac, 'afra.fra', 'units'  , 'm^2' )
575   CALL flioputa ( il_frac, 'afra.fra', 'title', 'lmdz-t mesh fraction, set to 1')
576
577   CALL flioputv ( il_frac, 'tlmd.fra', xfra )
578   CALL flioputv ( il_frac, 'aone.fra', 1.0_rl + 0.0_rl*xfra )
579   CALL flioputv ( il_frac, 'afra.fra', xfra )
580
581   CALL flioclo ( il_frac)
582   
583   !
584   IF ( l_read ) THEN 
585      !! Sortie fichier bathy
586      OPEN ( unit = 27, file = 'bathy.lmdz', FORM = 'formatted', STATUS = 'unknown', ACTION = 'write', POSITION = 'rewind' )
587      DO jaj = jpaj, 1, -1
588         WRITE ( unit = 27, FMT = '(92I1)' )   INT ( xmas(:,jaj), KIND = i_4)
589      END DO
590      CLOSE ( unit = 27)
591      !! Sortie fichier bathy flux iceberg
592      flux_iceberg = 0.0
593      DO jai = 1, jpai
594         DO jaj = 1, jpaj
595            IF ( xmas(jai,jaj) == 0 ) THEN
596               flux_iceberg(jai,jaj  ) = 2.0
597               flux_iceberg(jai,jaj+1) = 2.0
598               flux_iceberg(jai,jaj+2) = 1.0
599               flux_iceberg(jai,jaj+3) = 1.0
600               EXIT
601            END IF
602         END DO
603      END DO
604      !!
605      OPEN ( unit = 27, file = 'flux_iceberg', FORM = 'formatted', STATUS = 'unknown', ACTION = 'write', POSITION = 'rewind' )
606      DO jaj = jpaj, 1, -1
607         WRITE ( unit = 27, FMT = '(96F3.0)' )  flux_iceberg(:,jaj)
608      END DO
609   END IF
610   !!
611   STOP
612CONTAINS
613   ELEMENTAL FUNCTION clo_lon ( zlon, zlon0) ! Find closest longitude
614      IMPLICIT NONE
615      REAL (kind=rl) :: clo_lon
616      REAL (kind=rl), INTENT ( in) :: zlon, zlon0
617      REAL (kind=rl) :: zz, zp, zm, z0
618      INTEGER (kind=il) :: jn
619      z0 = zlon
620      DO jn = 1_il, 2_il
621         zz = ABS ( (z0           ) - zlon0 )
622         zp = ABS ( (z0 + 360.0_rl) - zlon0 )
623         zm = ABS ( (z0 - 360.0_rl) - zlon0 )
624         IF (      zp < MIN ( zm, zz) ) THEN
625            z0 = z0 + 360.0_rl
626         ELSE IF ( zm < MIN ( zp, zz) ) THEN
627            z0 = z0 - 360.0_rl
628         END IF
629      END DO
630      clo_lon = z0
631   END FUNCTION clo_lon
632
633   SUBROUTINE surface ( c_type )
634      IMPLICIT NONE
635      CHARACTER (len=*), OPTIONAL :: c_type
636      CHARACTER (len=8) :: cl_type
637      IF ( PRESENT (c_type)) THEN
638         cl_type = TRIM (c_type)
639      ELSE
640         cl_type = 'globe'
641      END IF
642
643      IF ( TRIM(c_type) == 'globe') THEN
644         WRITE (nout,*) 'Calcul des surfaces atm partout'
645         xare (:,:) = (xae(:,:)-xaw(:,:))*(yan(:,:)-yas(:,:)) * COS(rad*xlatt(:,:)) * rad * rad * ra * ra
646      END IF
647     
648      WRITE (nout,*) 'Calcul des surfaces atm aux poles'
649      xare (:, 1  ) = (xae(:,   1)-xaw(:,  1 ))*(1.0_rl - ABS(SIN(rad*yas(:,  1 )))) * rad * ra * ra
650      xare (:,jpaj) = (xae(:,jpaj)-xaw(:,jpaj))*(1.0_rl - ABS(SIN(rad*yan(:,jpaj)))) * rad * ra * ra
651
652   END SUBROUTINE surface
653
654   SUBROUTINE lmd1d_2d ( ztab1, ztab2 )
655      REAL (kind=rl), INTENT (in) , DIMENSION (:)   :: ztab1
656      REAL (kind=rl), INTENT (out), DIMENSION (:,:) :: ztab2
657      INTEGER :: jn, jn_min, jn_max
658
659      jn_min =  9999999
660      jn_max = -9999999
661
662      ztab2 ( :,  1  ) = ztab1 ( 1)
663      ztab2 ( :, jpai) = ztab1 ( nlmd)
664     
665      DO jaj = 2, jpai-1
666         DO jai = 1, jpai
667            jn = 1 + (jaj-2)*jpai + jai
668            jn_min = MIN ( jn, jn_min ) 
669            jn_max = MIN ( jn, jn_max ) 
670            ztab2 (jai, jaj) = ztab1 ( jn )
671         END DO
672      END DO
673      WRITE (nout,*) 'lmd1d_2d : ', jn_min, jn_max
674   END SUBROUTINE lmd1d_2d
675           
676END PROGRAM lmdgrille
Note: See TracBrowser for help on using the repository browser.