Changeset 280 for trunk/LMDZ.MARS/util

Timestamp:

Sep 6, 2011, 11:02:38 AM (13 years ago)

Author:

emillour

Message:

Update utilitaire "concatnc" du GCM Mars.
EM

File:

• trunk/LMDZ.MARS/util/concatnc.F90 (modified) (24 diffs)

trunk/LMDZ.MARS/util/concatnc.F90

@@ -10 +10 @@
 !   to output file (E. Millour, september 2006)
 !   if available (F. Forget, october 2006)
+! + handle 1D data (EM, January 2007)
 ! + ap(), bp()  (F. Forget, February 2008)
+! + handle the possibility that number of GCM layers (aps, bps
+!   or sigma) may be different from number of vertical levels
+!   of data (which is the case for outputs from 'zrecast')
+!   (EM, April 2010)
 ! ********************************************************
 
@@ -17 +22 @@
 include "netcdf.inc" ! NetCDF definitions
 
-character (len=50), dimension(200) :: file
+character (len=50), dimension(1000) :: file
 ! file(): input file(s) names(s)
 character (len=30), dimension(15) :: notconcat
@@ -43 +48 @@
 integer :: varid
 ! varid: [netcdf] variable ID #
-integer :: memolatlen,memolonlen,memoaltlen
+integer :: memolatlen=0,memolonlen=0,memoaltlen=0
 ! memolatlen: # of elements of lat(), read from the first input file
 ! memolonlen: # of elements of lon(), read from the first input file
@@ -62 +67 @@
 ! altdim: [netcdf] "altdim" dim ID
 ! timedim: [netcdf] "timedim" dim ID
+integer :: gcmlayerdim ! NetCDF dimension ID for # of layers in GCM
 integer :: latvar,lonvar,altvar,timevar
 ! latvar: [netcdf] ID of "latitude" variable
@@ -72 +78 @@
 ! altvar: # of elements of alt() array
 ! timelen: # of elemnets of time() array
-integer :: nout,latdimout,londimout,altdimout,apdimout,timedimout,timevarout
+integer :: GCM_layers ! number of GCM atmospheric layers (may not be
+! same as altlen if file is an output of zrecast)
+integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
 ! nout: [netcdf] output file ID
 ! latdimout: [netcdf] output latitude (dimension) ID
@@ -79 +87 @@
 ! timedimout: [netcdf] output time (dimension) ID
 ! timevarout: [netcdf] ID of output "Time" variable
+integer :: layerdimout ! NetCDF dimension ID for # of layers in GCM
+integer :: interlayerdimout ! dimension ID for # of interlayers in GCM
 integer :: reptime,rep,varidout
 ! reptime: total length of concatenated time() arrays
@@ -308 +318 @@
 !  write(*,*) "altlen: ",altlen
 
+! load size of aps() or sigma() (in case it is not altlen)
+   ! default is that GCM_layers=altlen
+   ! but for outputs of zrecast, it may be a different value
+   ierr=NF_INQ_DIMID(nid,"GCM_layers",gcmlayerdim)
+   if (ierr.ne.NF_NOERR) then
+     ! didn't find a GCM_layers dimension; therefore we have:
+     GCM_layers=altlen
+   else
+     ! load value of GCM_layers
+     ierr=NF_INQ_DIMLEN(nid,gcmlayerdim,GCM_layers)
+   endif
+!   write(*,*) "GCM_layers=",GCM_layers
+
 !==============================================================================
 ! 2.3. Read (and check compatibility of) dimensions of
@@ -330 +353 @@
 #endif
    ! Initialize output file's lat,lon,alt and time dimensions
-      call initiate (filename,lat,lon,alt,nout,&
-           latdimout,londimout,altdimout,apdimout,timedimout,timevarout)
+      call initiate (filename,lat,lon,alt,GCM_layers,nout,&
+       latdimout,londimout,altdimout,timedimout,&
+       layerdimout,interlayerdimout,timevarout)
    ! Initialize output file's aps,bps,ap,bp and phisinit variables
-     call init2(nid,lonlen,latlen,altlen,&
-                nout,londimout,latdimout,altdimout,apdimout)
+     call init2(nid,lonlen,latlen,altlen,GCM_layers,&
+                nout,londimout,latdimout,altdimout,&
+                layerdimout,interlayerdimout)
    
    else ! Not a first call,
@@ -424 +449 @@
       ! build dim(),corner() and edges() arrays
       ! and allocate var3d() array
-      if (ndim==3) then
+      if (ndim==1) then
+         allocate(var3d(timelen,1,1,1))
+         dim(1)=timedimout
+
+         ! start indexes (where data values will be written)
+         corner(1)=reptime+1
+         corner(2)=1
+         corner(3)=1
+         corner(4)=1
+
+         ! length (along dimensions) of block of data to be written
+         edges(1)=timelen 
+         edges(2)=1 
+         edges(3)=1
+         edges(4)=1
+
+      else if (ndim==3) then
          allocate(var3d(lonlen,latlen,timelen,1))
          dim(1)=londimout
@@ -527 +568 @@
 
 !******************************************************************************
-Subroutine initiate (filename,lat,lon,alt,&
-         nout,latdimout,londimout,altdimout,apdimout,timedimout,timevarout)
+subroutine initiate (filename,lat,lon,alt,GCM_layers,nout,&
+         latdimout,londimout,altdimout,timedimout,&
+         layerdimout,interlayerdimout,timevarout)
 !==============================================================================
 ! Purpose:
@@ -552 +594 @@
 real, dimension(:), intent(in):: alt
 ! alt(): altitude
+integer,intent(in) :: GCM_layers ! number of GCM layers
 integer, intent(out):: nout
 ! nout: [netcdf] file ID
@@ -560 +603 @@
 integer, intent(out):: altdimout
 ! altdimout: [netcdf] alt()  ID
-integer, intent(out):: apdimout
-! apdimout: [netcdf] ap()  ID
 integer, intent(out):: timedimout
 ! timedimout: [netcdf] "Time"  ID
+integer,intent(out) :: layerdimout
+! layerdimout: [netcdf] "GCM_layers" ID
+integer,intent(out) :: interlayerdimout
+! layerdimout: [netcdf] "GCM_layers+1" ID
 integer, intent(out):: timevarout
 ! timevarout: [netcdf] "Time" (considered as a variable) ID
@@ -593 +638 @@
 ierr = NF_DEF_DIM(nout, "longitude", size(lon), londimout)
 ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
-ierr = NF_DEF_DIM(nout, "interlayer",(size(alt)+1), apdimout)
 ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
+ierr = NF_DEF_DIM(nout, "GCM_layers", GCM_layers, layerdimout)
+ierr = NF_DEF_DIM(nout, "GCM_interlayers",GCM_layers+1,interlayerdimout)
 
 ! End netcdf define mode
@@ -660 +706 @@
 end Subroutine initiate
 !******************************************************************************
-subroutine init2(infid,lonlen,latlen,altlen, &
-                 outfid,londimout,latdimout,altdimout,apdimout)
+subroutine init2(infid,lonlen,latlen,altlen,GCM_layers, &
+                 outfid,londimout,latdimout,altdimout, &
+                 layerdimout,interlayerdimout)
 !==============================================================================
 ! Purpose:
-! Copy aps(), bps() and phisinit() from input file to outpout file
-! Copy ap(), bp() from input file to outpout file
+! Copy ap() , bp(), aps(), bps() and phisinit() from input file to outpout file
 !==============================================================================
 ! Remarks:
@@ -682 +728 @@
 integer, intent(in) :: latlen ! # of grid points along latitude
 integer, intent(in) :: altlen ! # of grid points along latitude
+integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
 integer, intent(in) :: outfid ! NetCDF output file ID
 integer, intent(in) :: londimout ! longitude dimension ID
 integer, intent(in) :: latdimout ! latitude dimension ID
 integer, intent(in) :: altdimout ! altitude dimension ID
-integer, intent(in) :: apdimout ! interlayer dimension ID
+integer, intent(in) :: layerdimout ! GCM_layers dimension ID
+integer, intent(in) :: interlayerdimout ! GCM_layers+1 dimension ID
 !==============================================================================
 ! Local variables:
@@ -692 +740 @@
 real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
 real,dimension(:),allocatable :: ap,bp ! hybrid vertical coordinates
+real,dimension(:),allocatable :: sigma ! sigma levels
 real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
-integer :: apsid,bpsid,phisinitid
+integer :: apsid,bpsid,sigmaid,phisinitid
 integer :: apid,bpid
 integer :: ierr
 integer :: tmpvarid ! temporary variable ID
 logical :: phis ! is "phisinit" available ?
-logical :: hybrid ! are "aps" "bps" "ap" "bp" available ?
+logical :: hybrid ! are "aps" and "bps" available ?
 
 !==============================================================================
@@ -705 +754 @@
 
 ! hybrid coordinate aps
-allocate(aps(altlen))
+!write(*,*) "aps: altlen=",altlen," GCM_layers=",GCM_layers
+allocate(aps(GCM_layers),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "init2: failed to allocate aps!"
+  stop
+endif
 ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
 if (ierr.ne.NF_NOERR) then
-  write(*,*) "Ooops. Failed to get aps ID. OK."
+  write(*,*) "Ooops. Failed to get aps ID. OK, will look for sigma coord."
   hybrid=.false.
 else
@@ -714 +768 @@
   hybrid=.true.
   if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed reading aps"
-  endif
-endif
-
-! hybrid coordinate bps
-allocate(bps(altlen))
-ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
+    stop "init2 Error: Failed reading aps"
+  endif
+
+  ! hybrid coordinate bps
+!  write(*,*) "bps: altlen=",altlen," GCM_layers=",GCM_layers
+  allocate(bps(GCM_layers),stat=ierr)
+  if (ierr.ne.0) then
+    write(*,*) "init2: failed to allocate bps!"
+    stop
+  endif
+  ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed to get bps ID."
+  endif
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed reading bps"
+  endif
+endif
+
+! hybrid coordinate ap
+allocate(ap(GCM_layers+1),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "init2: failed to allocate ap!"
+  stop
+else
+  ierr=NF_INQ_VARID(infid,"ap",tmpvarid)
+  if (ierr.ne.NF_NOERR) then
+    write(*,*) "Ooops. Failed to get ap ID. OK."
+    hybrid=.false.
+  else
+    ierr=NF_GET_VAR_REAL(infid,tmpvarid,ap)
+    hybrid=.true.
+    if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed reading ap"
+    endif
+  endif
+endif
+
+! hybrid coordinate bp
+allocate(bp(GCM_layers+1),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "init2: failed to allocate bp!"
+  stop
+else
+  ierr=NF_INQ_VARID(infid,"bp",tmpvarid)
+  if (ierr.ne.NF_NOERR) then
+    write(*,*) "Ooops. Failed to get bp ID. OK."
+    hybrid=.false.
+  else
+    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bp)
+    hybrid=.true.
+    if (ierr.ne.NF_NOERR) then
+      stop "Error: Failed reading bp"
+    endif
+  endif
+endif
+
+! sigma levels (if any)
+if (.not.hybrid) then
+  allocate(sigma(GCM_layers),stat=ierr)
+  if (ierr.ne.0) then
+    write(*,*) "init2: failed to allocate sigma"
+    stop
+  endif
+  ierr=NF_INQ_VARID(infid,"sigma",tmpvarid)
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed reading sigma"
+  endif
+endif ! of if (.not.hybrid)
+
+! ground geopotential phisinit
+allocate(phisinit(lonlen,latlen),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "init2: failed to allocate phisinit!"
+  stop
+endif
+ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
 if (ierr.ne.NF_NOERR) then
-  write(*,*) "Ooops. Failed to get bps ID. OK."
-  hybrid=.false.
-else
-  ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
-  hybrid=.true.
-  if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed reading bps"
-  endif
-endif
-
-! hybrid coordinate ap
-allocate(ap(altlen+1))
-ierr=NF_INQ_VARID(infid,"ap",tmpvarid)
-if (ierr.ne.NF_NOERR) then
-  write(*,*) "Ooops. Failed to get ap ID. OK."
-  hybrid=.false.
-else
-  ierr=NF_GET_VAR_REAL(infid,tmpvarid,ap)
-  hybrid=.true.
-  if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed reading ap"
-  endif
-endif
-
-! hybrid coordinate bp
-allocate(bp(altlen+1))
-ierr=NF_INQ_VARID(infid,"bp",tmpvarid)
-if (ierr.ne.NF_NOERR) then
-  write(*,*) "Ooops. Failed to get bp ID. OK."
-  hybrid=.false.
-else
-  ierr=NF_GET_VAR_REAL(infid,tmpvarid,bp)
-  hybrid=.true.
-  if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed reading bp"
-  endif
-endif
-
-! ground geopotential phisinit
-! ground geopotential phisinit
-ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
-allocate(phisinit(lonlen,latlen))
-if (ierr.ne.NF_NOERR) then
-  write(*,*) "Failed to get phisinit ID. We must be reading a ""stats"" file"
-  phisinit = 0.
+  write(*,*)"init2 warning: Failed to get phisinit ID."
   phis = .false.
 else
   ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
   if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed reading phisinit"
+    stop "init2 Error: Failed reading phisinit"
   endif
   phis = .true.
@@ -781 +863 @@
 
 !==============================================================================
-! 2.2. Hybrid coordinates aps() and bps()
+! 2.2. Hybrid coordinates ap() , bp(), aps() and bps()
 !==============================================================================
 if(hybrid) then 
 ! define aps
   call def_var(nout,"aps","hybrid pressure at midlayers"," ",1,&
-             (/altdimout/),apsid,ierr)
-  if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed to def_var aps"
+             (/layerdimout/),apsid,ierr)
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed to def_var aps"
   endif
 
@@ -798 +880 @@
 #endif
   if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed to write aps"
+    stop "init2 Error: Failed to write aps"
   endif
 
 ! define bps
   call def_var(nout,"bps","hybrid sigma at midlayers"," ",1,&
-             (/altdimout/),bpsid,ierr)
-  if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed to def_var bps"
+             (/layerdimout/),bpsid,ierr)
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed to def_var bps"
   endif
 
@@ -815 +897 @@
 #endif
   if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed to write bps"
+    stop "init2 Error: Failed to write bps"
   endif
 
 ! define ap
   call def_var(nout,"ap","hybrid sigma at interlayers"," ",1,&
-             (/apdimout/),apid,ierr)
+             (/interlayerdimout/),apid,ierr)
   if (ierr.ne.NF_NOERR) then
     stop "Error: Failed to def_var ap"
@@ -837 +919 @@
 ! define bp
   call def_var(nout,"bp","hybrid sigma at interlayers"," ",1,&
-             (/apdimout/),bpid,ierr)
+             (/interlayerdimout/),bpid,ierr)
   if (ierr.ne.NF_NOERR) then
     stop "Error: Failed to def_var bp"
@@ -851 +933 @@
     stop "Error: Failed to write bp"
   endif
-endif
+
+else
+! define sigma
+  call def_var(nout,"sigma","sigma at midlayers"," ",1,&
+             (/layerdimout/),sigmaid,ierr)
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed to def_var sigma"
+  endif
+! write sigma
+#ifdef NC_DOUBLE
+  ierr=NF_PUT_VAR_DOUBLE(outfid,sigmaid,sigma)
+#else
+  ierr=NF_PUT_VAR_REAL(outfid,sigmaid,sigma)
+#endif
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed to write sigma"
+  endif
+endif ! of if (hybrid)
 
 !==============================================================================
@@ -859 +958 @@
 IF (phis) THEN
 
-!define phisinit
- call def_var(nout,"phisinit","Ground level geopotential"," ",2,&
+  !define phisinit
+   call def_var(nout,"phisinit","Ground level geopotential"," ",2,&
             (/londimout,latdimout/),phisinitid,ierr)
- if (ierr.ne.NF_NOERR) then
-     stop "Error: Failed to def_var phisinit"
-  endif
-
-! write phisinit
-#ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(outfid,phisinitid,phisinit)
-#else
-ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
-#endif
-if (ierr.ne.NF_NOERR) then
-  stop "Error: Failed to write phisinit"
-endif
-
-ENDIF
+   if (ierr.ne.NF_NOERR) then
+     stop "init2 Error: Failed to def_var phisinit"
+   endif
+
+  ! write phisinit
+#ifdef NC_DOUBLE
+  ierr=NF_PUT_VAR_DOUBLE(outfid,phisinitid,phisinit)
+#else
+  ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
+#endif
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed to write phisinit"
+  endif
+
+ENDIF ! of IF (phis)
 
 
 ! Cleanup
-deallocate(aps)
-deallocate(bps)
-deallocate(ap)
-deallocate(bp)
-deallocate(phisinit)
+if (allocated(aps)) deallocate(aps)
+if (allocated(bps)) deallocate(bps)
+if (allocated(ap)) deallocate(ap)
+if (allocated(bp)) deallocate(bp)
+if (allocated(sigma)) deallocate(sigma)
+if (allocated(phisinit)) deallocate(phisinit)
 
 end subroutine init2