Changeset 290 for codes/icosagcm/branches

Timestamp:

10/22/14 10:21:38 (10 years ago)

Author:

ymipsl

Message:

Add vertical axis value to the restart file in order to be able to visualize the data file with VTK_mapper.

YM

File:

• codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/restart.f90 (modified) (5 diffs)

codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/restart.f90

@@ -21 +21 @@
   USE spherical_geom_mod
   USE transfert_mod
+  USE disvert_mod
  
   IMPLICIT NONE 
@@ -38 +39 @@
   CHARACTER(LEN=255) :: restart_file_name
   INTEGER,PARAMETER  :: nvert=6
-  INTEGER    ::  ncid, cellId, levId, edgeId,  vertid, lonId, latId, bounds_lonId, bounds_latId, nqId
+  INTEGER    ::  ncid, cellId, levId, edgeId,  vertid, lonId, latId, bounds_lonId, bounds_latId, nqId, levAxisId
   INTEGER    :: ind,ind_glo,i,j,k,nf
   INTEGER    :: status
@@ -91 +92 @@
       status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ vertId,cellId /),bounds_lonId)
       status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ vertId,cellId /),bounds_latId)
+      status = NF90_DEF_VAR(ncid,'lev',NF90_DOUBLE,(/ levId /),levAxisId)
+      status = NF90_PUT_ATT(ncid,levAxisId,"axis","Z")
+      status = NF90_PUT_ATT(ncid,levAxisId,"units","Pa")
+      status = NF90_PUT_ATT(ncid,levAxisId,"positive","down")
       
       DO nf=1,nfield
@@ -97 +102 @@
           IF (field(1)%ndim==2) THEN
             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(field(1)%name)),NF90_DOUBLE,(/ cellId /),fieldId(nf))
+            status = NF90_PUT_ATT(ncid,FieldId(nf),"coordinates","lon lat")
           ELSE IF (field(1)%ndim==3) THEN
             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(field(1)%name)),NF90_DOUBLE,(/ cellId, levId /),fieldId(nf))
+            status = NF90_PUT_ATT(ncid,FieldId(nf),"coordinates","lev lon lat")
           ELSE IF (field(1)%ndim==4) THEN
             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(field(1)%name)),NF90_DOUBLE,(/ cellId, levId,nqId /),fieldId(nf))
+            status = NF90_PUT_ATT(ncid,FieldId(nf),"coordinates","nq lev lon lat")
           ENDIF
-          status = NF90_PUT_ATT(ncid,FieldId(nf),"coordinates","lon lat")
         ELSE IF (field(1)%field_type==field_U) THEN
           IF (field(1)%ndim==2) THEN
@@ -139 +146 @@
       status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell_glo /))
       status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell_glo /))
+      status=NF90_PUT_VAR(ncid,levAxisId,REAL(presnivs,r8),start=(/ 1 /),count=(/ llm /))
     ENDIF