Index: trunk/LMDZ.MARS/util/concatnc.F90 =================================================================== --- trunk/LMDZ.MARS/util/concatnc.F90 (revision 275) +++ trunk/LMDZ.MARS/util/concatnc.F90 (revision 280) @@ -10,5 +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,5 +22,5 @@ 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,5 +48,5 @@ 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,4 +67,5 @@ ! 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,5 +78,7 @@ ! 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,4 +87,6 @@ ! 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,4 +318,17 @@ ! 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,9 +353,11 @@ #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,5 +449,21 @@ ! 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,6 +568,7 @@ !****************************************************************************** -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,4 +594,5 @@ 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,8 +603,10 @@ 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,6 +638,7 @@ 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,10 +706,10 @@ 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,9 +728,11 @@ 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,11 +740,12 @@ 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,8 +754,13 @@ ! 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,62 +768,90 @@ 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,12 +863,12 @@ !============================================================================== -! 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,12 +880,12 @@ #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,10 +897,10 @@ #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,5 +919,5 @@ ! 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,5 +933,22 @@ 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,30 +958,31 @@ 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