source: trunk/LMDZ.COMMON/libf/evolution/iostart_PEM.F90 @ 3296

Last change on this file since 3296 was 3206, checked in by jbclement, 10 months ago

PEM:
Cleanings of unused variables/arguments and bad type conversions.
JBC

File size: 35.7 KB
Line 
1MODULE iostart_PEM
2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3!!!
4!!! Purpose: Read and write specific netcdf for the PEM
5!!!
6!!!
7!!! Author: LL, inspired by iostart from the PCM
8!!!
9!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10
11    IMPLICIT NONE
12PRIVATE
13    INTEGER,SAVE :: nid_start ! NetCDF file identifier for startfi.nc file
14    INTEGER,SAVE :: nid_restart ! NetCDF file identifier for restartfi.nc file
15   
16    ! restartfi.nc file dimension identifiers: (see open_restartphy())
17    INTEGER,SAVE :: idim1 ! "index" dimension
18    INTEGER,SAVE :: idim2 ! "physical_points" dimension
19    INTEGER,SAVE :: idim3 ! "subsurface_layers" dimension
20    INTEGER,SAVE :: idim4 ! "nlayer_plus_1" dimension
21    INTEGER,SAVE :: idim5 ! "number_of_advected_fields" dimension
22    INTEGER,SAVE :: idim6 ! "nlayer" dimension
23    INTEGER,SAVE :: idim7 ! "Time" dimension
24    INTEGER,SAVE :: idim8 ! "nslope" dimension
25    INTEGER,SAVE :: idim9 ! "inter slope" dimension
26    INTEGER,SAVE :: timeindex ! current time index (for time-dependent fields)
27    INTEGER,PARAMETER :: length=100 ! size of tab_cntrl array
28   
29    INTERFACE get_field
30      MODULE PROCEDURE Get_field_r1,Get_field_r2,Get_field_r3
31    END INTERFACE get_field
32   
33    INTERFACE get_var
34      MODULE PROCEDURE get_var_r0,Get_var_r1,Get_var_r2,Get_var_r3
35    END INTERFACE get_var
36
37    INTERFACE put_field
38      MODULE PROCEDURE put_field_r1,put_field_r2,put_field_r3
39    END INTERFACE put_field
40
41    INTERFACE put_var
42      MODULE PROCEDURE put_var_r0,put_var_r1,put_var_r2,put_var_r3
43    END INTERFACE put_var
44
45    PUBLIC nid_start, length
46    PUBLIC get_field,get_var,put_field,put_var
47    PUBLIC inquire_dimension, inquire_dimension_length
48    PUBLIC inquire_field, inquire_field_ndims
49    PUBLIC open_startphy,close_startphy,open_restartphy,close_restartphy
50   
51CONTAINS
52
53  SUBROUTINE open_startphy(filename)
54  USE netcdf, only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE, nf90_strerror
55  USE mod_phys_lmdz_para, only: is_master, bcast
56  IMPLICIT NONE
57    CHARACTER(LEN=*) :: filename
58    INTEGER          :: ierr
59
60    IF (is_master) THEN
61      ierr = NF90_OPEN (filename, NF90_NOWRITE, nid_start)
62      IF (ierr.NE.NF90_NOERR) THEN
63        write(*,*)'open_startphy: problem opening file '//trim(filename)
64        write(*,*)trim(nf90_strerror(ierr))
65        CALL abort_physic("open_startphy","Cannot open file",1)
66      ENDIF
67    ENDIF
68   
69    CALL bcast(nid_start) ! tell all procs about nid_start
70 
71  END SUBROUTINE open_startphy
72
73  SUBROUTINE close_startphy
74  USE netcdf, only: NF90_CLOSE
75  USE mod_phys_lmdz_para, only: is_master
76  IMPLICIT NONE
77    INTEGER          :: ierr
78
79    IF (is_master) THEN
80        ierr = NF90_CLOSE (nid_start)
81    ENDIF
82
83  END SUBROUTINE close_startphy
84
85
86  FUNCTION inquire_field(Field_name)
87  ! check if a given field is present in the input file
88  USE netcdf, only: NF90_INQ_VARID, NF90_NOERR
89  USE mod_phys_lmdz_para, only: is_master, bcast
90  IMPLICIT NONE
91    CHARACTER(LEN=*),INTENT(IN) :: Field_name
92    LOGICAL :: inquire_field
93    INTEGER :: varid
94    INTEGER :: ierr
95   
96    IF (is_master) THEN
97      ierr=NF90_INQ_VARID(nid_start,Field_name,varid)
98      IF (ierr==NF90_NOERR) THEN
99        Inquire_field=.TRUE.
100      ELSE
101        Inquire_field=.FALSE.
102      ENDIF
103    ENDIF
104
105    CALL bcast(inquire_field)
106
107  END FUNCTION inquire_field
108
109
110  FUNCTION inquire_field_ndims(Field_name)
111  ! give the number of dimensions of "Field_name" stored in the input file
112  USE netcdf, only: nf90_inq_varid, nf90_inquire_variable, &
113                    NF90_NOERR, nf90_strerror
114  USE mod_phys_lmdz_para, only: is_master, bcast
115  IMPLICIT NONE
116    CHARACTER(LEN=*),INTENT(IN) :: Field_name
117    INTEGER :: inquire_field_ndims
118    INTEGER :: varid
119    INTEGER :: ierr
120   
121    IF (is_master) THEN
122      ierr=nf90_inq_varid(nid_start,Field_name,varid)
123      ierr=nf90_inquire_variable(nid_start,varid,&
124                                  ndims=inquire_field_ndims)
125      IF (ierr.NE.NF90_NOERR) THEN
126        write(*,*)'inquire_field_ndims: problem obtaining ndims of '&
127                  //trim(field_name)
128        write(*,*)trim(nf90_strerror(ierr))
129        CALL abort_physic("inquire_field_ndims","Failed to get ndims",1)
130      ENDIF
131    ENDIF
132
133    CALL bcast(inquire_field_ndims)
134
135  END FUNCTION inquire_field_ndims
136
137
138  FUNCTION inquire_dimension(Field_name)
139  ! check if a given dimension is present in the input file
140  USE netcdf, only: nf90_inq_dimid, NF90_NOERR
141  USE mod_phys_lmdz_para, only: is_master, bcast
142  IMPLICIT NONE
143    CHARACTER(LEN=*),INTENT(IN) :: Field_name
144    LOGICAL :: inquire_dimension
145    INTEGER :: varid
146    INTEGER :: ierr
147   
148    IF (is_master) THEN
149      ierr=NF90_INQ_DIMID(nid_start,Field_name,varid)
150      IF (ierr==NF90_NOERR) THEN
151        Inquire_dimension=.TRUE.
152      ELSE
153        Inquire_dimension=.FALSE.
154      ENDIF
155    ENDIF
156
157    CALL bcast(inquire_dimension)
158
159  END FUNCTION inquire_dimension
160
161  FUNCTION inquire_dimension_length(Field_name)
162  ! give the length of the "Field_name" dimension stored in the input file
163  USE netcdf, only: nf90_inquire_dimension, nf90_inq_dimid, &
164                    NF90_NOERR, nf90_strerror
165  USE mod_phys_lmdz_para, only: is_master, bcast
166  IMPLICIT NONE
167    CHARACTER(LEN=*),INTENT(IN) :: Field_name
168    INTEGER :: inquire_dimension_length
169    INTEGER :: varid
170    INTEGER :: ierr
171   
172    IF (is_master) THEN
173      ierr=nf90_inq_dimid(nid_start,Field_name,varid)
174      ierr=nf90_inquire_dimension(nid_start,varid,&
175                                  len=inquire_dimension_length)
176      IF (ierr.NE.NF90_NOERR) THEN
177        write(*,*)'inquire_field_length: problem obtaining length of '&
178                  //trim(field_name)
179        write(*,*)trim(nf90_strerror(ierr))
180        CALL abort_physic("inquire_field_ndims","Failed to get length",1)
181      ENDIF
182    ENDIF
183
184    CALL bcast(inquire_dimension_length)
185
186  END FUNCTION inquire_dimension_length
187
188
189
190  SUBROUTINE Get_Field_r1(field_name,field,found,timeindex)
191  ! For a surface field
192  use mod_grid_phy_lmdz, only: klon_glo ! number of atmospheric columns (full grid)
193  IMPLICIT NONE
194    CHARACTER(LEN=*),INTENT(IN)    :: Field_name
195    REAL,INTENT(INOUT)               :: Field(:)
196    LOGICAL,INTENT(OUT),OPTIONAL   :: found
197    INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
198
199    integer :: edges(4), corners(4)
200
201    edges(1)=klon_glo
202    edges(2:4)=1
203    corners(1:4)=1
204    if (PRESENT(timeindex)) then
205      corners(2)=timeindex
206    endif
207
208    IF (PRESENT(found)) THEN
209      CALL Get_field_rgen(field_name,field,1,corners,edges,found)
210    ELSE
211      CALL Get_field_rgen(field_name,field,1,corners,edges)
212    ENDIF
213     
214  END SUBROUTINE Get_Field_r1
215 
216  SUBROUTINE Get_Field_r2(field_name,field,found,timeindex)
217  ! For a "3D" horizontal-vertical field
218  use mod_grid_phy_lmdz, only: klon_glo ! number of atmospheric columns (full grid)
219  IMPLICIT NONE
220    CHARACTER(LEN=*),INTENT(IN)    :: Field_name
221    REAL,INTENT(INOUT)               :: Field(:,:)
222    LOGICAL,INTENT(OUT),OPTIONAL   :: found
223    INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
224
225    integer :: edges(4), corners(4)
226
227    edges(1)=klon_glo
228    edges(2)=size(field,2)
229    edges(3:4)=1
230    corners(1:4)=1
231    if (PRESENT(timeindex)) then
232      corners(3)=timeindex
233    endif
234   
235    IF (PRESENT(found)) THEN
236      CALL Get_field_rgen(field_name,field,size(field,2),&
237                          corners,edges,found)
238    ELSE
239      CALL Get_field_rgen(field_name,field,size(field,2),&
240                          corners,edges)
241    ENDIF
242
243     
244  END SUBROUTINE Get_Field_r2
245 
246  SUBROUTINE Get_Field_r3(field_name,field,found,timeindex)
247  ! for a "4D" field surf/alt/??
248  use mod_grid_phy_lmdz, only: klon_glo ! number of atmospheric columns (full grid)
249  IMPLICIT NONE
250    CHARACTER(LEN=*),INTENT(IN)    :: Field_name
251    REAL,INTENT(INOUT)               :: Field(:,:,:)
252    LOGICAL,INTENT(OUT),OPTIONAL   :: found
253    INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
254
255    integer :: edges(4), corners(4)
256
257    edges(1)=klon_glo
258    edges(2)=size(field,2)
259    edges(3)=size(field,3)
260    edges(4)=1
261    corners(1:4)=1
262    if (PRESENT(timeindex)) then
263      corners(4)=timeindex
264    endif
265   
266    IF (PRESENT(found)) THEN
267      CALL Get_field_rgen(field_name,field,size(field,2)*size(field,3),&
268                          corners,edges,found)
269    ELSE
270      CALL Get_field_rgen(field_name,field,size(field,2)*size(field,3),&
271                          corners,edges)
272    ENDIF
273     
274  END SUBROUTINE Get_Field_r3
275 
276  SUBROUTINE Get_field_rgen(field_name,field,field_size, &
277                            corners,edges,found)
278  USE netcdf
279  USE dimphy
280  USE mod_grid_phy_lmdz
281  USE mod_phys_lmdz_para
282  IMPLICIT NONE
283    CHARACTER(LEN=*) :: Field_name
284    INTEGER          :: field_size
285    REAL             :: field(klon,field_size)
286    INTEGER,INTENT(IN) :: corners(4)
287    INTEGER,INTENT(IN) :: edges(4)
288    LOGICAL,OPTIONAL :: found
289   
290    REAL    :: field_glo(klon_glo,field_size)
291    LOGICAL :: tmp_found
292    INTEGER :: varid
293    INTEGER :: ierr
294   
295    IF (is_master) THEN
296 
297      ierr=NF90_INQ_VARID(nid_start,Field_name,varid)
298     
299      IF (ierr==NF90_NOERR) THEN
300        CALL body(field_glo)
301        tmp_found=.TRUE.
302      ELSE
303        tmp_found=.FALSE.
304      ENDIF
305   
306    ENDIF
307   
308    CALL bcast(tmp_found)
309
310    IF (tmp_found) THEN
311      CALL scatter(field_glo,field)
312    ENDIF
313   
314    IF (PRESENT(found)) THEN
315      found=tmp_found
316    ELSE
317      IF (.NOT. tmp_found) THEN
318        write(*,*) 'get_field_rgen: Field <'//field_name//'> not found'
319        CALL abort_physic("get_field_rgen","Failed to get field",1)
320      ENDIF
321    ENDIF
322 
323   
324    CONTAINS
325     
326     SUBROUTINE body(field_glo)
327       REAL :: field_glo(klon_glo*field_size)
328         ierr=NF90_GET_VAR(nid_start,varid,field_glo,corners,edges)
329         IF (ierr/=NF90_NOERR) THEN
330           ! La variable exist dans le fichier mais la lecture a echouee.
331           write(*,*) 'get_field_rgen: Failed reading <'//field_name//'>'
332
333!           IF (field_name=='CLWCON' .OR. field_name=='RNEBCON' .OR. field_name=='RATQS') THEN
334!              ! Essaye de lire le variable sur surface uniqument, comme fait avant
335!              field_glo(:)=0.
336!              ierr=NF90_GET_VAR(nid_start,varid,field_glo(1:klon_glo))
337!              IF (ierr/=NF90_NOERR) THEN
338!                 write(*,*) 'phyetat0: Lecture echouee aussi en 2D pour <'//field_name//'>'
339!                 CALL abort
340!              ELSE
341!                 write(*,*) 'phyetat0: La variable <'//field_name//'> lu sur surface seulement'!, selon ancien format, le reste mis a zero'
342!              END IF
343!           ELSE
344              CALL abort_physic("get_field_rgen","Failed to read field",1)
345!           ENDIF
346         ENDIF
347
348     END SUBROUTINE body
349
350  END SUBROUTINE Get_field_rgen
351
352
353  SUBROUTINE get_var_r0(var_name,var,found)
354  ! Get a scalar from input file
355  IMPLICIT NONE 
356    CHARACTER(LEN=*),INTENT(IN)  :: var_name
357    REAL,INTENT(INOUT)             :: var
358    LOGICAL,OPTIONAL,INTENT(OUT) :: found
359
360    REAL                         :: varout(1)
361   
362    IF (PRESENT(found)) THEN
363      CALL Get_var_rgen(var_name,varout,size(varout),found)
364    ELSE
365      CALL Get_var_rgen(var_name,varout,size(varout))
366    ENDIF
367    var=varout(1)
368 
369  END SUBROUTINE get_var_r0
370
371  SUBROUTINE get_var_r1(var_name,var,found)
372  ! Get a vector from input file
373  IMPLICIT NONE 
374    CHARACTER(LEN=*),INTENT(IN)  :: var_name
375    REAL,INTENT(INOUT)             :: var(:)
376    LOGICAL,OPTIONAL,INTENT(OUT) :: found
377   
378    IF (PRESENT(found)) THEN
379      CALL Get_var_rgen(var_name,var,size(var),found)
380    ELSE
381      CALL Get_var_rgen(var_name,var,size(var))
382    ENDIF
383 
384  END SUBROUTINE get_var_r1
385
386  SUBROUTINE get_var_r2(var_name,var,found)
387  ! Get a 2D field from input file
388  IMPLICIT NONE 
389    CHARACTER(LEN=*),INTENT(IN)  :: var_name
390    REAL,INTENT(OUT)             :: var(:,:)
391    LOGICAL,OPTIONAL,INTENT(OUT) :: found
392   
393    IF (PRESENT(found)) THEN
394      CALL Get_var_rgen(var_name,var,size(var),found)
395    ELSE
396      CALL Get_var_rgen(var_name,var,size(var))
397    ENDIF
398 
399  END SUBROUTINE get_var_r2
400
401  SUBROUTINE get_var_r3(var_name,var,found)
402  ! Get a 3D field frominput file
403  IMPLICIT NONE 
404    CHARACTER(LEN=*),INTENT(IN)  :: var_name
405    REAL,INTENT(INOUT)             :: var(:,:,:)
406    LOGICAL,OPTIONAL,INTENT(OUT) :: found
407   
408    IF (PRESENT(found)) THEN
409      CALL Get_var_rgen(var_name,var,size(var),found)
410    ELSE
411      CALL Get_var_rgen(var_name,var,size(var))
412    ENDIF
413 
414  END SUBROUTINE get_var_r3
415
416  SUBROUTINE Get_var_rgen(var_name,var,var_size,found)
417  USE netcdf
418  USE dimphy
419  USE mod_grid_phy_lmdz
420  USE mod_phys_lmdz_para
421  IMPLICIT NONE
422    CHARACTER(LEN=*) :: var_name
423    INTEGER          :: var_size
424    REAL             :: var(var_size)
425    LOGICAL,OPTIONAL :: found
426   
427    LOGICAL :: tmp_found
428    INTEGER :: varid
429    INTEGER :: ierr
430   
431    IF (is_mpi_root .AND. is_omp_root) THEN
432 
433      ierr=NF90_INQ_VARID(nid_start,var_name,varid)
434     
435      IF (ierr==NF90_NOERR) THEN
436        ierr=NF90_GET_VAR(nid_start,varid,var)
437        IF (ierr/=NF90_NOERR) THEN
438          write(*,*) 'phyetat0: Failed loading <'//trim(var_name)//'>'
439          CALL abort_physic("get_var_rgen","Failed to read variable",1)
440        ENDIF
441        tmp_found=.TRUE.
442      ELSE
443        tmp_found=.FALSE.
444      ENDIF
445   
446    ENDIF
447   
448    CALL bcast(tmp_found)
449
450    IF (tmp_found) THEN
451      CALL bcast(var)
452    ENDIF
453   
454    IF (PRESENT(found)) THEN
455      found=tmp_found
456    ELSE
457      IF (.NOT. tmp_found) THEN
458        write(*,*) 'phyetat0: Variable <'//trim(var_name)//'> not found'
459        CALL abort_physic("get_var_rgen","Failed to read variable",1)
460      ENDIF
461    ENDIF
462
463  END SUBROUTINE Get_var_rgen
464
465!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
466
467  SUBROUTINE open_restartphy(filename)
468  USE netcdf, only: NF90_CREATE, NF90_CLOBBER, NF90_64BIT_OFFSET, &
469                    NF90_NOERR, nf90_strerror, &
470                    NF90_PUT_ATT, NF90_GLOBAL, NF90_DEF_DIM, &
471                    NF90_UNLIMITED, NF90_ENDDEF, &
472                    NF90_WRITE, NF90_OPEN
473  USE mod_phys_lmdz_para, only: is_master
474  USE mod_grid_phy_lmdz, only: klon_glo
475  USE dimphy, only: klev, klevp1
476#ifndef CPP_STD
477  USE tracer_mod, only: nqmx
478#else
479  use tracer_h, only:   nqtot
480#endif
481  USE comsoil_h_PEM, only:  nsoilmx_PEM
482  USE comslope_mod, only: nslope
483  IMPLICIT NONE
484    CHARACTER(LEN=*),INTENT(IN) :: filename
485    INTEGER                     :: ierr
486    LOGICAL,SAVE :: already_created=.false.
487
488#ifdef CPP_STD
489    INTEGER :: nqmx
490    nqmx=nqtot
491#endif
492   
493    IF (is_master) THEN
494      if (.not.already_created) then
495        ! At the very first call, create the file
496        ierr=NF90_CREATE(filename,IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), &
497                          nid_restart)
498        IF (ierr/=NF90_NOERR) THEN
499          write(*,*)'open_restartphy: problem creating file '//trim(filename)
500          write(*,*)trim(nf90_strerror(ierr))
501          CALL abort_physic("open_restartphy","Failed creating file",1)
502        ENDIF
503        already_created=.true.
504      else
505        ! Just open the file
506        ierr=NF90_OPEN(filename,NF90_WRITE,nid_restart)
507        IF (ierr/=NF90_NOERR) THEN
508          write(*,*)'open_restartphy: problem opening file '//trim(filename)
509          write(*,*)trim(nf90_strerror(ierr))
510          CALL abort_physic("open_restartphy","Failed opening file",1)
511        ENDIF
512        return
513      endif ! of if (.not.already_created)
514
515      ierr=NF90_PUT_ATT(nid_restart,NF90_GLOBAL,"title",&
516                        "Physics start file")
517      IF (ierr/=NF90_NOERR) THEN
518        write(*,*)'open_restartphy: problem writing title '
519        write(*,*)trim(nf90_strerror(ierr))
520      ENDIF
521
522      ierr=NF90_DEF_DIM(nid_restart,"index",length,idim1)
523      IF (ierr/=NF90_NOERR) THEN
524        write(*,*)'open_restartphy: problem defining index dimension '
525        write(*,*)trim(nf90_strerror(ierr))
526        CALL abort_physic("open_restartphy","Failed defining index",1)
527      ENDIF
528     
529      ierr=NF90_DEF_DIM(nid_restart,"physical_points",klon_glo,idim2)
530      IF (ierr/=NF90_NOERR) THEN
531        write(*,*)'open_restartphy: problem defining physical_points dimension '
532        write(*,*)trim(nf90_strerror(ierr))
533        CALL abort_physic("open_restartphy","Failed defining physical_points",1)
534      ENDIF
535     
536      ierr=NF90_DEF_DIM(nid_restart,"subsurface_layers",nsoilmx_PEM,idim3)
537      IF (ierr/=NF90_NOERR) THEN
538        write(*,*)'open_restartphy: problem defining subsurface_layers dimension '
539        write(*,*)trim(nf90_strerror(ierr))
540        CALL abort_physic("open_restartphy","Failed defining subsurface_layers",1)
541      ENDIF
542     
543      ierr=NF90_DEF_DIM(nid_restart,"nlayer_plus_1",klevp1,idim4)
544      IF (ierr/=NF90_NOERR) THEN
545        write(*,*)'open_restartphy: problem defining nlayer_plus_1 dimension '
546        write(*,*)trim(nf90_strerror(ierr))
547        CALL abort_physic("open_restartphy","Failed defining nlayer_plus_1",1)
548      ENDIF
549     
550      if (nqmx.ne.0) then
551        ierr=NF90_DEF_DIM(nid_restart,"number_of_advected_fields",nqmx,idim5)
552      else
553        ! pretend nqmx=1 because 0 size implies unlimited dimension for NetCDF
554        ierr=NF90_DEF_DIM(nid_restart,"number_of_advected_fields",1,idim5)
555      endif
556      IF (ierr/=NF90_NOERR) THEN
557        write(*,*)'open_restartphy: problem defining number_of_advected_fields dimension '
558        write(*,*)trim(nf90_strerror(ierr))
559        CALL abort_physic("open_restartphy","Failed defining number_of_advected_fields",1)
560      ENDIF
561
562      ierr=NF90_DEF_DIM(nid_restart,"nlayer",klev,idim6)
563      IF (ierr/=NF90_NOERR) THEN
564        write(*,*)'open_restartphy: problem defining nlayer dimension '
565        write(*,*)trim(nf90_strerror(ierr))
566        CALL abort_physic("open_restartphy","Failed defining nlayer",1)
567      ENDIF
568     
569      ierr=NF90_DEF_DIM(nid_restart,"Time",NF90_UNLIMITED,idim7)
570      IF (ierr/=NF90_NOERR) THEN
571        write(*,*)'open_restartphy: problem defining Time dimension '
572        write(*,*)trim(nf90_strerror(ierr))
573        CALL abort_physic("open_restartphy","Failed defining Time",1)
574      ENDIF
575
576      ierr=NF90_DEF_DIM(nid_restart,"nslope",nslope,idim8)
577      IF (ierr/=NF90_NOERR) THEN
578        write(*,*)'phyredem: problem defining nslope dimension'
579        write(*,*)trim(nf90_strerror(ierr))
580        CALL ABORT
581      ENDIF
582
583      ierr=NF90_DEF_DIM(nid_restart,"inter slope",nslope+1,idim9)
584      IF (ierr/=NF90_NOERR) THEN
585        write(*,*)'phyredem: problem defining inter slope dimension'
586        write(*,*)trim(nf90_strerror(ierr))
587        CALL ABORT
588      ENDIF
589
590
591      ierr=NF90_ENDDEF(nid_restart)
592      IF (ierr/=NF90_NOERR) THEN
593        write(*,*)'open_restartphy: problem ending definition mode '
594        write(*,*)trim(nf90_strerror(ierr))
595        CALL abort_physic("open_restartphy","Failed ending definition mode",1)
596      ENDIF
597    ENDIF
598
599  END SUBROUTINE open_restartphy
600
601  SUBROUTINE close_restartphy
602  USE netcdf, only: NF90_CLOSE
603  USE mod_phys_lmdz_para, only: is_master
604  IMPLICIT NONE
605    INTEGER          :: ierr
606
607    IF (is_master) THEN
608      ierr = NF90_CLOSE (nid_restart)
609    ENDIF
610 
611  END SUBROUTINE close_restartphy
612
613  SUBROUTINE put_field_r1(field_name,title,field,time)
614  ! For a surface field
615  IMPLICIT NONE
616  CHARACTER(LEN=*),INTENT(IN)    :: field_name
617  CHARACTER(LEN=*),INTENT(IN)    :: title
618  REAL,INTENT(IN)                :: field(:)
619  REAL,OPTIONAL,INTENT(IN)       :: time
620 
621  IF (present(time)) THEN
622    ! if timeindex is present, it is a time-dependent variable
623    CALL put_field_rgen(field_name,title,field,1,time)
624  ELSE
625    CALL put_field_rgen(field_name,title,field,1)
626  ENDIF
627 
628  END SUBROUTINE put_field_r1
629
630  SUBROUTINE put_field_r2(field_name,title,field,time)
631  ! For a "3D" horizontal-vertical field
632  IMPLICIT NONE
633  CHARACTER(LEN=*),INTENT(IN)    :: field_name
634  CHARACTER(LEN=*),INTENT(IN)    :: title
635  REAL,INTENT(IN)                :: field(:,:)
636  REAL,OPTIONAL,INTENT(IN)       :: time
637 
638  IF (present(time)) THEN
639    ! if timeindex is present, it is a time-dependent variable
640    CALL put_field_rgen(field_name,title,field,size(field,2),time)
641  ELSE
642    CALL put_field_rgen(field_name,title,field,size(field,2))
643  ENDIF
644 
645  END SUBROUTINE put_field_r2
646
647  SUBROUTINE put_field_r3(field_name,title,field,time)
648  ! For a "4D" field surf/alt/??
649  IMPLICIT NONE
650  CHARACTER(LEN=*),INTENT(IN)    :: field_name
651  CHARACTER(LEN=*),INTENT(IN)    :: title
652  REAL,INTENT(IN)                :: field(:,:,:)
653  REAL,OPTIONAL,INTENT(IN)       :: time
654 
655  IF (present(time)) THEN
656    ! if timeindex is present, it is a time-dependent variable
657    CALL put_field_rgen(field_name,title,field,size(field,2)*size(field,3),&
658                        time)
659  ELSE 
660    CALL put_field_rgen(field_name,title,field,size(field,2)*size(field,3))
661  ENDIF
662 
663  END SUBROUTINE put_field_r3
664 
665  SUBROUTINE put_field_rgen(field_name,title,field,field_size,time)
666  USE netcdf
667  USE dimphy
668  USE comsoil_h_PEM, only:  nsoilmx_PEM
669  USE comslope_mod, ONLY: nslope
670  USE mod_grid_phy_lmdz
671  USE mod_phys_lmdz_para
672  IMPLICIT NONE
673  CHARACTER(LEN=*),INTENT(IN)    :: field_name
674  CHARACTER(LEN=*),INTENT(IN)    :: title
675  INTEGER,INTENT(IN)             :: field_size
676  REAL,INTENT(IN)                :: field(klon,field_size)
677  REAL,OPTIONAL,INTENT(IN)       :: time
678 
679  REAL                           :: field_glo(klon_glo,field_size)
680  INTEGER                        :: ierr
681  INTEGER                        :: nvarid
682   
683    CALL gather(field,field_glo)
684   
685    IF (is_master) THEN
686
687      IF (field_size==1) THEN
688        ! input is a 1D "surface field" array
689        if (.not.present(time)) then ! for a time-independent field
690          ierr=NF90_REDEF(nid_restart)
691#ifdef NC_DOUBLE
692          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
693                            (/idim2/),nvarid)
694#else
695          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
696                            (/idim2/),nvarid)
697#endif
698          if (ierr.ne.NF90_NOERR) then
699            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
700            write(*,*)trim(nf90_strerror(ierr))
701          endif
702          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
703          ierr=NF90_ENDDEF(nid_restart)
704          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo)
705        else
706          ! check if the variable has already been defined:
707          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
708          if (ierr/=NF90_NOERR) then ! variable not found, define it
709            ierr=NF90_REDEF(nid_restart)
710#ifdef NC_DOUBLE
711            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
712                              (/idim2,idim7/),nvarid)
713#else
714            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
715                              (/idim2,idim7/),nvarid)
716#endif
717            if (ierr.ne.NF90_NOERR) then
718              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
719              write(*,*)trim(nf90_strerror(ierr))
720            endif
721            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
722            ierr=NF90_ENDDEF(nid_restart)
723          endif
724          ! Write the variable
725          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
726                            start=(/1,timeindex/))
727        endif ! of if (.not.present(timeindex))
728
729      ELSE IF (field_size==klev) THEN
730        ! input is a 2D "atmospheric field" array
731        if (.not.present(time)) then ! for a time-independent field
732          ierr=NF90_REDEF(nid_restart)
733#ifdef NC_DOUBLE
734          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
735                            (/idim2,idim6/),nvarid)
736#else
737          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
738                            (/idim2,idim6/),nvarid)
739#endif
740          if (ierr.ne.NF90_NOERR) then
741            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
742            write(*,*)trim(nf90_strerror(ierr))
743          endif
744          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
745          ierr=NF90_ENDDEF(nid_restart)
746          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo)
747        else
748          ! check if the variable has already been defined:
749          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
750          if (ierr/=NF90_NOERR) then ! variable not found, define it
751            ierr=NF90_REDEF(nid_restart)
752#ifdef NC_DOUBLE
753            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
754                              (/idim2,idim6,idim7/),nvarid)
755#else
756            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
757                              (/idim2,idim6,idim7/),nvarid)
758#endif
759            if (ierr.ne.NF90_NOERR) then
760              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
761              write(*,*)trim(nf90_strerror(ierr))
762            endif
763            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
764            ierr=NF90_ENDDEF(nid_restart)
765          endif
766          ! Write the variable
767          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
768                            start=(/1,1,timeindex/))
769        endif ! of if (.not.present(time))
770
771      ELSE IF (field_size==klevp1) THEN
772        ! input is a 2D "interlayer atmospheric field" array
773        if (.not.present(time)) then ! for a time-independent field
774          ierr=NF90_REDEF(nid_restart)
775#ifdef NC_DOUBLE
776          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
777                            (/idim2,idim4/),nvarid)
778#else
779          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
780                            (/idim2,idim4/),nvarid)
781#endif
782          if (ierr.ne.NF90_NOERR) then
783            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
784            write(*,*)trim(nf90_strerror(ierr))
785          endif
786          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
787          ierr = NF90_ENDDEF(nid_restart)
788          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
789        else
790          ! check if the variable has already been defined:
791          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
792          if (ierr/=NF90_NOERR) then ! variable not found, define it
793            ierr=NF90_REDEF(nid_restart)
794#ifdef NC_DOUBLE
795            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
796                              (/idim2,idim4,idim7/),nvarid)
797#else
798            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
799                              (/idim2,idim4,idim7/),nvarid)
800#endif
801            if (ierr.ne.NF90_NOERR) then
802              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
803              write(*,*)trim(nf90_strerror(ierr))
804            endif
805            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
806            ierr=NF90_ENDDEF(nid_restart)
807          endif
808          ! Write the variable
809          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
810                            start=(/1,1,timeindex/))
811        endif ! of if (.not.present(timeindex))
812
813      ELSE IF (field_size==nsoilmx_PEM) THEN
814        ! input is a 2D "subsurface field" array
815        if (.not.present(time)) then ! for a time-independent field
816          ierr = NF90_REDEF(nid_restart)
817#ifdef NC_DOUBLE
818          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
819                            (/idim2,idim3/),nvarid)
820#else
821          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
822                            (/idim2,idim3/),nvarid)
823#endif
824          if (ierr.ne.NF90_NOERR) then
825            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
826            write(*,*)trim(nf90_strerror(ierr))
827          endif
828          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
829          ierr = NF90_ENDDEF(nid_restart)
830          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
831        else
832          ! check if the variable has already been defined:
833          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
834          if (ierr/=NF90_NOERR) then ! variable not found, define it
835            ierr=NF90_REDEF(nid_restart)
836#ifdef NC_DOUBLE
837            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
838                              (/idim2,idim3,idim7/),nvarid)
839#else
840            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
841                              (/idim2,idim3,idim7/),nvarid)
842#endif
843           if (ierr.ne.NF90_NOERR) then
844              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
845              write(*,*)trim(nf90_strerror(ierr))
846            endif
847            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
848            ierr=NF90_ENDDEF(nid_restart)
849          endif
850          ! Write the variable
851          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
852                            start=(/1,1,timeindex/))
853
854        endif ! of if (.not.present(time))
855
856
857      ELSE IF (field_size==nslope) THEN
858        ! input is a 2D "subsurface field" array
859        if (.not.present(time)) then ! for a time-independent field
860          ierr = NF90_REDEF(nid_restart)
861#ifdef NC_DOUBLE
862          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
863                            (/idim2,idim8/),nvarid)
864#else
865          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
866                            (/idim2,idim8/),nvarid)
867#endif
868          if (ierr.ne.NF90_NOERR) then
869            write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
870            write(*,*)trim(nf90_strerror(ierr))
871          endif
872          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
873          ierr = NF90_ENDDEF(nid_restart)
874          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
875        else
876          ! check if the variable has already been defined:
877          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
878          if (ierr/=NF90_NOERR) then ! variable not found, define it
879            ierr=NF90_REDEF(nid_restart)
880#ifdef NC_DOUBLE
881            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
882                              (/idim2,idim8,idim7/),nvarid)
883#else
884            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
885                              (/idim2,idim8,idim7/),nvarid)
886#endif
887           if (ierr.ne.NF90_NOERR) then
888              write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
889              write(*,*)trim(nf90_strerror(ierr))
890            endif
891            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
892            ierr=NF90_ENDDEF(nid_restart)
893          endif
894          ! Write the variable
895          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
896                            start=(/1,1,timeindex/))
897
898        endif ! of if (.not.present(time))
899
900      ELSE
901        write(*,*) "Error phyredem(put_field_rgen) : wrong dimension for ",trim(field_name)
902        write(*,*) "  field_size =",field_size
903        CALL abort_physic("put_field_rgen","wrong field dimension",1)
904      ENDIF
905
906      ! Check the writting of field to file went OK
907      if (ierr.ne.NF90_NOERR) then
908        write(*,*) " Error phyredem(put_field_rgen) : failed writing ",trim(field_name)
909        write(*,*)trim(nf90_strerror(ierr))
910        CALL abort_physic("put_field_rgen","Failed writing field",1)
911      endif
912
913    ENDIF ! of IF (is_master)
914   
915  END SUBROUTINE put_field_rgen 
916 
917  SUBROUTINE put_var_r0(var_name,title,var)
918  ! Put a scalar in file
919   IMPLICIT NONE
920     CHARACTER(LEN=*),INTENT(IN) :: var_name
921     CHARACTER(LEN=*),INTENT(IN) :: title
922     REAL,INTENT(IN)             :: var
923     REAL                        :: varin(1)
924     
925     varin(1)=var
926     
927     CALL put_var_rgen(var_name,title,varin,size(varin))
928
929  END SUBROUTINE put_var_r0
930
931
932  SUBROUTINE put_var_r1(var_name,title,var)
933  ! Put a vector in file
934   IMPLICIT NONE
935     CHARACTER(LEN=*),INTENT(IN) :: var_name
936     CHARACTER(LEN=*),INTENT(IN) :: title
937     REAL,INTENT(IN)             :: var(:)
938     
939     CALL put_var_rgen(var_name,title,var,size(var))
940
941  END SUBROUTINE put_var_r1
942 
943  SUBROUTINE put_var_r2(var_name,title,var)
944  ! Put a 2D field in file
945   IMPLICIT NONE
946     CHARACTER(LEN=*),INTENT(IN) :: var_name
947     CHARACTER(LEN=*),INTENT(IN) :: title
948     REAL,INTENT(IN)             :: var(:,:)
949     
950     CALL put_var_rgen(var_name,title,var,size(var))
951
952  END SUBROUTINE put_var_r2     
953 
954  SUBROUTINE put_var_r3(var_name,title,var)
955  ! Put a 3D field in file
956   IMPLICIT NONE
957     CHARACTER(LEN=*),INTENT(IN) :: var_name
958     CHARACTER(LEN=*),INTENT(IN) :: title
959     REAL,INTENT(IN)             :: var(:,:,:)
960     
961     CALL put_var_rgen(var_name,title,var,size(var))
962
963  END SUBROUTINE put_var_r3
964
965  SUBROUTINE put_var_rgen(var_name,title,var,var_size)
966  USE netcdf, only: NF90_REDEF, NF90_DEF_VAR, NF90_ENDDEF, NF90_PUT_VAR, &
967                    NF90_FLOAT, NF90_DOUBLE, &
968                    NF90_PUT_ATT, NF90_NOERR, nf90_strerror, &
969                    nf90_inq_dimid, nf90_inquire_dimension, NF90_INQ_VARID
970  USE comsoil_h_PEM, only:  nsoilmx_PEM
971  USE comslope_mod, only: nslope
972  USE mod_phys_lmdz_para, only: is_master
973  IMPLICIT NONE
974     CHARACTER(LEN=*),INTENT(IN) :: var_name
975     CHARACTER(LEN=*),INTENT(IN) :: title
976     INTEGER,INTENT(IN)          :: var_size
977     REAL,INTENT(IN)             :: var(var_size)
978     
979     INTEGER :: ierr
980     INTEGER :: nvarid
981     INTEGER :: idim1d
982     logical,save :: firsttime=.true.
983         
984    IF (is_master) THEN
985
986      IF (var_name=="Time") THEN
987        ! Very specific case of "Time" variable
988        if (firsttime) then
989          ! Create the "Time variable"
990          ierr=NF90_REDEF(nid_restart)
991#ifdef NC_DOUBLE
992          ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_DOUBLE,&
993                            (/idim7/),nvarid)
994#else
995          ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_FLOAT,&
996                            (/idim7/),nvarid)
997#endif
998          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
999          ierr=NF90_ENDDEF(nid_restart)
1000         
1001          firsttime=.false.
1002        endif
1003        ! Append "time" value
1004        ! get current length of "Time" dimension
1005        ierr=nf90_inq_dimid(nid_restart,var_name,idim1d)
1006        ierr=nf90_inquire_dimension(nid_restart,idim1d,len=timeindex)
1007        timeindex=timeindex+1
1008        ierr=NF90_INQ_VARID(nid_restart,var_name,nvarid)
1009        ierr=NF90_PUT_VAR(nid_restart,nvarid,var,&
1010                           start=(/timeindex/))
1011        IF (ierr/=NF90_NOERR) THEN
1012          write(*,*)'put_var_rgen: problem writing Time'
1013          write(*,*)trim(nf90_strerror(ierr))
1014          CALL abort_physic("get_var_rgen","Failed to write Time",1)
1015        ENDIF
1016        return ! nothing left to do
1017      ELSEIF (var_size==length) THEN
1018        ! We know it is a "controle" kind of 1D array
1019        idim1d=idim1
1020      ELSEIF (var_size==nsoilmx_PEM) THEN
1021        ! We know it is an  "mlayer" kind of 1D array
1022        idim1d=idim3
1023      ELSEIF (var_size==nslope+1) THEN
1024        ! We know it is an  "inter slope" kind of 1D array
1025        idim1d=idim9
1026      ELSE
1027        write(*,*) "put_var_rgen error : wrong dimension"
1028        write(*,*) "  var_size =",var_size
1029        CALL abort_physic("get_var_rgen","Wrong variable dimension",1)
1030
1031      ENDIF ! of IF (var_size==length) THEN
1032
1033      ! Swich to NetCDF define mode
1034      ierr=NF90_REDEF (nid_restart)
1035      ! Define the variable
1036#ifdef NC_DOUBLE
1037      ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_DOUBLE,(/idim1d/),nvarid)
1038#else
1039      ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_FLOAT,(/idim1d/),nvarid)
1040#endif
1041      ! Add a "title" attribute
1042      IF (LEN_TRIM(title)>0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
1043      ! Swich out of define mode
1044      ierr=NF90_ENDDEF(nid_restart)
1045      ! Write variable to file
1046      ierr=NF90_PUT_VAR(nid_restart,nvarid,var)
1047      IF (ierr/=NF90_NOERR) THEN
1048        write(*,*)'put_var_rgen: problem writing '//trim(var_name)
1049        write(*,*)trim(nf90_strerror(ierr))
1050        CALL abort_physic("get_var_rgen","Failed writing variable",1)
1051      ENDIF
1052    ENDIF ! of IF (is_master)
1053   
1054  END SUBROUTINE put_var_rgen     
1055
1056END MODULE iostart_PEM
Note: See TracBrowser for help on using the repository browser.