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

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

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