Index: /BOL/Class_Reg/Makefile.linux =================================================================== --- /BOL/Class_Reg/Makefile.linux (revision 1585) +++ /BOL/Class_Reg/Makefile.linux (revision 1585) @@ -0,0 +1,5 @@ +LIBNETCDF=/u/fairhead/IPCC/CMOR/netcdf-3.5.1/lib +INCNETCDF=/u/fairhead/IPCC/CMOR/netcdf-3.5.1/include + +geo2reg : geo2reg.F90 + pgf90 -I${INCNETCDF} -L${LIBNETCDF} -lnetcdf geo2reg.F90 -o geo2reg Index: /BOL/Class_Reg/config.def =================================================================== --- /BOL/Class_Reg/config.def (revision 1585) +++ /BOL/Class_Reg/config.def (revision 1585) @@ -0,0 +1,8 @@ +# type de classification (w500 pour l'instant) +bintype=w500 +# sur quelle surface? (glob, ocea, mixt, land) +cnat=land +# limites des latitudes a traiter -lat0 < donnees < lat0 +lat0=30.0 +# unite de sortie +lunout = 6 Index: /BOL/Class_Reg/geo2reg.F90 =================================================================== --- /BOL/Class_Reg/geo2reg.F90 (revision 1585) +++ /BOL/Class_Reg/geo2reg.F90 (revision 1585) @@ -0,0 +1,700 @@ +! +! $Header$ +! +!------------------------------------------------------------------- +! Make relationship between any variable and w500 at resolution 2.5x2.5 +! -> mean seasonal cycle of the relationship over the period 1985-1989 +! -> interannual anomalies (by removing the mean seasonal cycle) +! -> decomposition of interannual anomalies into thermo, dynam and mix components +! +! Specify: +! - undefined value of the data +! - min/max/step of the w500 bins +! - nb of months in timeseries +! - montly index of January 1985 +! +! Sandrine Bony, March 2004 +! Repris pour faire du netcdf, L. Fairhead 2005/02 +!------------------------------------------------------------------- + + implicit none + +#include "netcdf.inc" + + integer imax,jmax,an,mmax,nbinmax,noc + parameter (imax=360,jmax=180,mmax=1700,an=12,nbinmax=200,noc=5) + + integer nrun + parameter (nrun=1) + + real zmsk(imax,jmax) + real surfmax,zsurf(imax,jmax) + logical ok_msk(imax,jmax) + + character*20 namegcm, run + + real surftest + real lat0, undef_x +! parameter (lat0 = 30.0) ! limits tropical belt + parameter (undef_x = 1.e+20 ) ! undefined value + + real, allocatable, dimension(:,:) :: surf, msk, w, x + real, allocatable, dimension(:,:) :: x_binw + real, allocatable, dimension(:,:) :: w_binw + real, allocatable, dimension(:,:) :: xpdf,nx_binw,sx_binw + real, allocatable, dimension(:) :: xpdfmean + + integer nb_bin, nbr_inf, nbr_sup + real min_bin, max_bin, step_bin, w_mult, x_mult + real xpdftot(mmax) + real xpdfmeantot + + + + INTEGER i,j,m,ir,im,jm,lm,itime,n, il, l + real undef, undef_w, pi, x1, w1, msk1 + integer :: iostat + character (len=512) :: line_read + character (len=20) :: first_part + character (len=128) :: second_part + character cnat*4, bintype*4 + +! netcdf stuff + real :: time + real, dimension(:), allocatable :: var_dim + INTEGER :: ierr, nvdims, lunout, in_file_id, ierrs, ndims, natts + integer :: lonid, latid, in_time_id, varid, out_file_id, ngatts, levid + integer :: nbin_id, var_bin_id, time_id, var_time_id, in_var_time_id + INTEGER :: var_pdf_id, var_nx_id, var_w_id, var_x_id, outlev_id + integer :: var_lev_id + integer, dimension(4) :: vdims, start, count + character (len=80) :: long_name, varname, attname + character (len=132) :: in_file + real, dimension(:), allocatable ::lat + + logical :: var_3d = .false. + + integer :: iargc + external iargc +!--------------------------------------------------------------------- +!-- General specifications : +!--------------------------------------------------------------------- + +!******************************* FIN INTERFACE ********************* +! Lecture du fichier de configuration + open (20, IOSTAT=iostat, file='config.def',form='formatted') + if (iostat /= 0) THEN + write(*,*)'Erreur ouverture du fichier config.def' + stop + endif + + config: do + read(20,'(A)',iostat=iostat)line_read + if (iostat /= 0) exit + line_read = trim(line_read) + IF (INDEX(line_read, '#') /= 1) THEN + first_part = trim(line_read(1:INDEX(line_read, '=')-1)) + second_part = trim(line_read(INDEX(line_read, '=')+1:)) + selectcase(first_part) + case('bintype') + bintype = trim(second_part) + case('cnat') + cnat = trim(second_part) + case('lunout') + read(second_part,'(i)') lunout + case('lat0') + read(second_part,'(f)') lat0 + end select + endif + enddo config + if (iostat > 0) then + write(lunout,*) & + & 'Probleme de lecture du fichier config.def, iostat = ',iostat + stop + endif + close(20) + +! undefined value for the output: + undef = 999.999 ! for output only + +!------------------------------------------------------------ +! Definition of the w500 bins: +!------------------------------------------------------------ + + if (bintype.eq.'w500') then + min_bin = -200.0 + max_bin = 200.0 + step_bin = 10.0 + step_bin = 5.0 + w_mult = 864.0 ! Pa/s -> hPa/day +! w_mult = 1.0 ! Pa/s -> hPa/day + endif + +! scale factor for surface temperature, precip and variable x: + x_mult = 1.0 + + undef_w = undef + +!--------------------------------------------------------------------- +!-- preliminaries: +!--------------------------------------------------------------------- + + if (cnat.ne.'ocea'.and.cnat.ne.'land'.and.cnat.ne.'glob'.and.cnat.ne.'mixt')& +& then + write(*,*) 'erreur cnat ',cnat + stop + endif + + nb_bin = (max_bin-min_bin)/step_bin + 1 + if (nb_bin.gt.nbinmax) then + write(*,*) 'augmenter le nb de bins' + write(*,*) 'nbinmax, nb_bin = ',nbinmax,nb_bin + stop + endif + +! Il faut deux arguments a l appel: fichier et variable + CALL getarg(1, in_file) + if (iargc() == 0 .OR. iargc() /= 2 ) THEN + write(lunout,*)' ' + write(lunout,*)' Utilisation de ce programme: ' + write(lunout,*)' ./geo2reg nom_de_fichier [variable]' + write(lunout,*) & + & ' ou nom_de_fichier est le nom du fichier a traiter' + write(lunout,*) & + & ' et variable la variable a traiter [optionel]' + write(lunout,*)' ' + write(lunout,*)' ./geo2reg -h sort ce message' + write(lunout,*)' ' + stop + endif + CALL getarg(2, varname) + +! Ouverture du fichier a traiter (histmth) + ierr = NF_OPEN(in_file, NF_NOwrite, in_file_id) + if (ierr /= NF_NOERR) then + write(lunout,*)NF_STRERROR(ierr) + stop + endif +! +! lire im, jm, lm et temps + ierrs = 0 + ierr = NF_INQ_DIMID(in_file_id, 'lon', lonid) + ierrs = ierrs + ierr + ierr = NF_INQ_DIMLEN(in_file_id, lonid, im) + ierrs = ierrs + ierr + if (ierrs /= 2 * NF_NOERR) THEN + write(lunout,*)'Pb. avec la lecture de la longitude' + stop + endif + ierrs = 0 + ierr = NF_INQ_DIMID(in_file_id, 'lat', latid) + ierrs = ierrs + ierr + ierr = NF_INQ_DIMLEN(in_file_id, latid, jm) + ierrs = ierrs + ierr + allocate (lat(jm)) + ierr = NF_INQ_VARID(in_file_id,'lat', latid) + ierrs = ierrs + ierr + ierr = NF_GET_VAR_REAL(in_file_id, latid, lat) + ierrs = ierrs + ierr + if (ierrs /= 4 * NF_NOERR) THEN + write(lunout,*)'Pb. avec la lecture de la latitude' + stop + endif + ierrs = 0 + ierr = NF_INQ_DIMID(in_file_id, 'presnivs', levid) + ierrs = ierrs + ierr + ierr = NF_INQ_DIMLEN(in_file_id, levid, lm) + ierrs = ierrs + ierr + if (ierrs /= 2 * NF_NOERR) THEN + write(lunout,*)'Pb. avec la lecture des niveaux verticaux' + stop + endif + + + ierrs = 0 + ierr = NF_INQ_DIMID(in_file_id, 'time_counter', in_time_id) + ierrs = ierrs + ierr + ierr = NF_INQ_DIMLEN(in_file_id, in_time_id, itime) + ierrs = ierrs + ierr + if (ierrs /= 2 * NF_NOERR) THEN + write(lunout,*)'Pb. avec la lecture du temps' + stop + endif + +! lecture de l aire des mailles + allocate (surf(im,jm)) + ierrs = 0 + ierr = NF_INQ_VARID(in_file_id, 'aire', varid) + ierrs = ierrs + ierr + ierr = NF_GET_VAR_REAL(in_file_id, varid, surf) + ierrs = ierrs + ierr + if (ierrs /= 2 * NF_NOERR) THEN + write(lunout,*)'Pb. avec la lecture de l''aire' + stop + endif + +!-- compute gridbox areas: + + + pi = acos(-1.) + surfmax=0. + surftest=0.0 + do i = 1, im + do j = 1, jm + surftest=surftest+surf(i,j) + if(surf(i,j).gt.surfmax) surfmax=surf(i,j) + enddo + enddo + write(*,*) 'sfce totale = 4pi ',surftest/pi + +!--------------------------------------------------------------- +! Lecture du masque +! + allocate (msk(im,jm)) + if (cnat.ne.'glob') then + ierrs = 0 + ierr = NF_INQ_VARID(in_file_id, 'pourc_ter', varid) + ierrs = ierrs + ierr + ierr = NF_GET_VAR_REAL(in_file_id, varid, msk) + ierrs = ierrs + ierr + if (ierrs /= 2 * NF_NOERR) THEN + write(lunout,*)'Pb. avec la lecture du masque' + stop + endif + else + do i = 1, im + do j = 1, jm + msk(i,j) = undef + enddo + enddo + endif + + + do i = 1, im + do j = 1, jm + ok_msk(i,j) = .FALSE. + if (cnat.eq.'ocea' .and. msk(i,j).le.0.01) ok_msk(i,j) = .TRUE. + if (cnat.eq.'mixt' .and. msk(i,j).gt.0.01.and. msk(i,j).le.0.99) & + & ok_msk(i,j) = .TRUE. + if (cnat.eq.'land' .and. msk(i,j).gt.0.99) ok_msk(i,j) = .TRUE. + if (cnat.eq.'glob') ok_msk(i,j) = .TRUE. + enddo + enddo + +! +! Pour savoir si la variable a traiter est 3D +! + ierrs = 0 + ierr = NF_INQ_VARID(in_file_id, varname, varid) + ierrs = ierrs + ierr + ierrs = NF_INQ_VARNDIMS(in_file_id, varid, ndims) + if (ierrs /= 2 * NF_NOERR) THEN + write(lunout,*)'Pb. avec la lecture de ', varname + stop + endif + if (ndims == 4) var_3d = .true. + +!-------------------------------------------------------------------- +! creation du fichier de sortie et entete +!-------------------------------------------------------------------- + + + ierr = NF_CREATE('out.nc', NF_CLOBBER, out_file_id) + if (ierr /= NF_NOERR) then + write(lunout,*)NF_STRERROR(ierr) + stop + endif + ierr = NF_INQ_NATTS(in_file_id, ngatts) + do i = 1, ngatts + ierr = NF_INQ_ATTNAME(in_file_id, NF_GLOBAL, i, attname) + ierr = NF_COPY_ATT(in_file_id, NF_GLOBAL, attname, out_file_id, NF_GLOBAL) + enddo +! +! definition des dimensions (nbin, level, temps) +! nbin: + ierrs = 0 + ierr = NF_DEF_DIM(out_file_id, 'nbin', nb_bin, nbin_id) + ierrs = ierrs + ierr + nvdims = 1 + vdims(1) = nbin_id + ierr = NF_DEF_VAR(out_file_id, 'nbin', NF_FLOAT, nvdims, vdims, var_bin_id) + ierrs = ierrs + ierr + ierr = NF_PUT_ATT_TEXT(out_file_id, var_bin_id, 'units', 5, 'hPa/j') + ierrs = ierrs + ierr + ierr = NF_PUT_ATT_TEXT(out_file_id, var_bin_id, 'long_name', 15, & + & 'number of bins') + ierrs = ierrs + ierr + if (ierrs /= 4 * NF_NOERR) THEN + write(lunout,*)'Pb. dans la definition de nbin' + stop + endif +! +! y bidon pour grads + ierr = NF_DEF_DIM(out_file_id, 'lat', 1, latid) + nvdims = 1 + vdims(1) = latid + ierr = NF_DEF_VAR(out_file_id, 'lat', NF_FLOAT, nvdims, vdims,latid) + ierr = NF_PUT_ATT_TEXT(out_file_id, latid, 'units', 13, 'degrees_north') + +! +! presnivs + ierrs = 0 + ierr = NF_DEF_DIM(out_file_id, 'presnivs', lm, outlev_id) + ierrs = ierrs + ierr + nvdims = 1 + vdims(1) = outlev_id + ierr = NF_DEF_VAR(out_file_id,'presnivs',NF_FLOAT, nvdims, vdims, var_lev_id) + ierrs = ierrs + ierr + ierr = NF_INQ_VARID(in_file_id, 'presnivs', levid) + ierrs = ierrs + ierr + ierr = NF_INQ_VARNATTS(in_file_id, levid, natts) + ierrs = ierrs + ierr + do i = 1, natts + attname = '' + ierr = NF_INQ_ATTNAME(in_file_id, levid, i, attname) + ierrs = ierrs + ierr + ierr = NF_COPY_ATT(in_file_id, levid, attname, out_file_id, var_lev_id) + ierrs = ierrs + ierr + enddo + IF (ierrs /= (4 + natts) * NF_NOERR) THEN + write(lunout,*)'Pb. dans la definition de presnivs' + stop + endif + +! +! temps: + ierrs = 0 + ierr = NF_DEF_DIM(out_file_id, 'time_counter', NF_UNLIMITED, time_id) + ierrs = ierrs + ierr + nvdims = 1 + vdims(1) = time_id + ierr = NF_DEF_VAR(out_file_id, 'time_counter', NF_FLOAT, nvdims, vdims, & + & var_time_id) + ierrs = ierrs + ierr + ierr = NF_INQ_VARID(in_file_id, 'time_counter', in_var_time_id) + if (ierr /= NF_NOERR) & + & ierr = NF_INQ_VARID(in_file_id, 't_ave_02592000', in_var_time_id) + if (ierr == NF_NOERR) then + ierr = NF_COPY_ATT(in_file_id, in_var_time_id, 'units', out_file_id, & + & var_time_id) + else + ierr = NF_PUT_ATT_TEXT(out_file_id, var_time_id, 'units', 34, & + & 'seconds since 1860-01-01 00:00:00') + in_var_time_id = 0 + endif + ierrs = ierrs + ierr + ierr = NF_PUT_ATT_TEXT(out_file_id, var_time_id, 'long_name', 9, & + & 'Time axis') + ierrs = ierrs + ierr + if (ierrs /= 4 * NF_NOERR) THEN + write(lunout,*)'Pb. dans la definition de time_counter' + stop + endif +! +! Definition des variables a ecrire: +! pdf: + ierrs = 0 + nvdims = 2 + vdims(1) = nbin_id + vdims(2) = time_id + ierr = NF_DEF_VAR(out_file_id, 'pdf', NF_FLOAT, nvdims, vdims, & + & var_pdf_id) + ierrs = ierrs + ierr + ierr = NF_PUT_ATT_TEXT(out_file_id, var_pdf_id, 'long_name', 16, & + & 'monthly PDF w500') + ierrs = ierrs + ierr + ierr = NF_PUT_ATT_REAL(out_file_id, var_pdf_id, 'missing_value', NF_FLOAT, & + & 1, undef) + ierrs = ierrs + ierr + if (ierrs /= 3 * NF_NOERR) THEN + write(lunout,*)'Pb. dans la definition de pdf' + stop + endif +! nx: + ierrs = 0 + nvdims = 2 + vdims(1) = nbin_id + vdims(2) = time_id + ierr = NF_DEF_VAR(out_file_id, 'nx', NF_FLOAT, nvdims, vdims, & + & var_nx_id) + ierrs = ierrs + ierr + ierr = NF_PUT_ATT_TEXT(out_file_id, var_nx_id, 'long_name', 16, & + & 'nb of points in w500 bin') + ierrs = ierrs + ierr + if (ierrs /= 2 * NF_NOERR) THEN + write(lunout,*)'Pb. dans la definition de nx' + stop + endif +! +! w: + ierrs = 0 + nvdims = 2 + vdims(1) = nbin_id + vdims(2) = time_id + ierr = NF_DEF_VAR(out_file_id, 'w', NF_FLOAT, nvdims, vdims, & + & var_w_id) + ierrs = ierrs + ierr + long_name = 'mean relationship w(binw)' + ierr = NF_PUT_ATT_TEXT(out_file_id, var_w_id, 'long_name', 25, & + & long_name) + ierrs = ierrs + ierr + ierr = NF_PUT_ATT_REAL(out_file_id, var_w_id, 'missing_value', NF_FLOAT, & + & 1, undef) + ierrs = ierrs + ierr + if (ierrs /= 3 * NF_NOERR) THEN + write(lunout,*)'Pb. dans la definition de w' + stop + endif +! +! x: + ierrs = 0 + nvdims = 2 + vdims(1) = nbin_id + vdims(2) = time_id + if (var_3d) then + nvdims = 3 + vdims(2) = outlev_id + vdims(3) = time_id + endif + ierr = NF_DEF_VAR(out_file_id, varname, NF_FLOAT, nvdims, vdims, & + & var_x_id) + ierrs = ierrs + ierr + long_name = 'mean relationship '//trim(varname) + ierr = NF_PUT_ATT_TEXT(out_file_id, var_x_id, 'long_name', len(long_name), & + & long_name) + ierrs = ierrs + ierr + ierr = NF_PUT_ATT_REAL(out_file_id, var_x_id, 'missing_value', NF_FLOAT, & + & 1, undef) + ierrs = ierrs + ierr + if (ierrs /= 3 * NF_NOERR) THEN + write(lunout,*)'Pb. dans la definition de ',varname + stop + endif + ierr = NF_ENDDEF(out_file_id) + +! +! On passe aux calculs +! +! initialisations: + + allocate(xpdfmean(nb_bin)) + allocate(xpdf(nb_bin, itime)) + allocate(nx_binw(nb_bin, itime)) + allocate(sx_binw(nb_bin, itime)) + allocate(w_binw(nb_bin, itime)) + allocate(x_binw(nb_bin, itime)) + +!-- temporal loop: + il = 1 + if (var_3d) then + il = lm + endif + + do l = 1, il + xpdfmeantot = 0. + xpdfmean = 0.0 + x_binw = 0.0 + xpdf = 0.0 + nx_binw = 0.0 + sx_binw = 0.0 + w_binw = 0.0 + xpdftot = 0. + do m = 1, itime ! loop over timeperiod of x_data + +! read data: + start(1) = 1 + start(2) = 1 + start(3) = m + count(1) = im + count(2) = jm + count(3) = 1 +! vit. verticale à 500hP + if (.not. allocated(w)) allocate (w(im,jm)) + ierrs = 0 + ierr = NF_INQ_VARID(in_file_id, 'w500', varid) + ierrs = ierrs + ierr + ierr = NF_GET_VARA_REAL(in_file_id, varid, start, count, w) + ierrs = ierrs + ierr + if (ierrs /= 2* NF_NOERR) THEN + write(lunout,*)'Pb. avec la lecture de la vit. vert. à 500hP' + stop + endif +! +! variable a traiter + if (.not. allocated(x)) allocate (x(im,jm)) + if (var_3d) then + start(3) = l + count(3) = 1 + start(4) = m + count(4) = 1 + endif + ierrs = 0 + ierr = NF_INQ_VARID(in_file_id, varname, varid) + ierrs = ierrs + ierr + ierr = NF_GET_VARA_REAL(in_file_id, varid, start, count, x) + ierrs = ierrs + ierr + if (ierrs /= 2 * NF_NOERR) THEN + write(lunout,*)'Pb. avec la lecture de ', varname + stop + endif +! + do j=1,jm + do i=1,im + zmsk(i,j)=0. + enddo + enddo + nbr_inf = 0; nbr_sup = 0 +! tropical belt: + do j = 1, jm + if (ABS(lat(j)).le.lat0) then +! loop over grid points: + do i = 1, im + if (x(i,j).ne.undef_x .and. w(i,j).ne.undef_w & + & .and. ok_msk(i,j) ) then + zmsk(i,j)=1. + w1 = w(i,j)*w_mult + x1 = x(i,j)*x_mult +! bin w500: + if (w1 < min_bin) then + ir = 1 + nbr_inf = nbr_inf + 1 + else if (w1 > max_bin) then + ir = nb_bin + nbr_sup = nbr_sup + 1 + else + ir = INT((w1-min_bin)/step_bin) + 1 + endif +! monthly PDF : + xpdfmeantot = xpdfmeantot + surf(i,j) + xpdfmean(ir) = xpdfmean(ir) + surf(i,j) + xpdftot(m) = xpdftot(m) + surf(i,j) + xpdf(ir,m) = xpdf(ir,m) + surf(i,j) +! monthly x-w relationship: + nx_binw(ir,m) = nx_binw(ir,m) + 1.0 + sx_binw(ir,m) = sx_binw(ir,m) + surf(i,j) + x_binw(ir,m) = x_binw(ir,m) + x1*surf(i,j) + w_binw(ir,m) = w_binw(ir,m) + w1*surf(i,j) + + endif + enddo ! i + endif ! lat + enddo ! j + if (nbr_inf /=0) & + & write(lunout,*)'nbre de points ou w500 < ',min_bin,' = ', nbr_inf + if (nbr_sup /=0) & + & write(lunout,*)'nbre de points ou w500 > ',max_bin,' = ', nbr_sup + enddo ! m + +! normalize PDF: + + do ir = 1, nb_bin + xpdfmean(ir) = xpdfmean(ir)/xpdfmeantot + enddo + + do m = 1, itime + if (xpdftot(m).gt.0.) then + do ir = 1, nb_bin + xpdf(ir,m) = xpdf(ir,m)/xpdftot(m) + enddo + else + do ir = 1, nb_bin + xpdf(ir,m) = undef + enddo + endif + enddo + + + do ir = 1, nb_bin + do m = 1, itime + if (nx_binw(ir,m).gt.1.) then +! if (nx_binw(ir,m).gt.0.) then + x_binw(ir,m) = x_binw(ir,m)/sx_binw(ir,m) + w_binw(ir,m) = w_binw(ir,m)/sx_binw(ir,m) + else + x_binw(ir,m) = undef + w_binw(ir,m) = undef + endif + enddo + enddo ! ir + +! +! ecriture du fichier sortie +! + + allocate(var_dim(nb_bin)) + do ir = 1, nb_bin + var_dim(ir) = min_bin + step_bin * (ir -1) + enddo + ierr = NF_PUT_VAR_REAL(out_file_id, var_bin_id, var_dim) + if (ierr /= NF_NOERR) then + write(lunout,*)NF_STRERROR(ierr) + stop + endif + deallocate(var_dim) +! +! ecriture y bidon +! + ierr = NF_PUT_VAR1_REAL(out_file_id, latid, 1, 0.) + +! +! ecriture niveaux verticaux +! + allocate(var_dim(lm)) + ierrs = 0 + ierr = NF_GET_VAR_REAL(in_file_id, levid, var_dim) + ierrs = ierrs + ierr + ierr = NF_PUT_VAR_REAL(out_file_id, var_lev_id, var_dim) + ierrs = ierrs + ierr + if (ierrs /= 2 * NF_NOERR) THEN + write(lunout,*)'Pb. d''ecriture niveaux verticaux' + stop + endif + deallocate(var_dim) + + + + do m = 1, itime + start(1) = 1 + start(2) = m + count(1) = nb_bin + count(2) = 1 + ierrs = 0 +! +! ecriture temps: + if (in_var_time_id == 0) then + time = (m - 1) * 86400. * 30. + else + ierr = NF_GET_VAR1_REAL(in_file_id, in_var_time_id, m, time) + endif + ierr = NF_PUT_VAR1_REAL(out_file_id, var_time_id, m, time) + ierrs = ierrs + ierr +! +! pdf, nb_bin, x, w: + ierr = NF_PUT_VARA_REAL(out_file_id, var_pdf_id, start, count, & + & xpdf(1,m)) + ierrs = ierrs + ierr + ierr = NF_PUT_VARA_REAL(out_file_id, var_nx_id, start, count, & + & nx_binw(1,m)) + ierrs = ierrs + ierr + ierr = NF_PUT_VARA_REAL(out_file_id, var_w_id, start, count, & + & w_binw(1,m)) + ierrs = ierrs + ierr + if (var_3d) then + start(2) = l + count(2) = 1 + start(3) = m + count(3) = 1 + endif + ierr = NF_PUT_VARA_REAL(out_file_id, var_x_id, start, count, & + & x_binw(1,m)) + ierrs = ierrs + ierr + if (ierrs /= 5 * NF_NOERR) THEN + write(lunout,*)'Pb. d''ecriture' + stop + endif + enddo ! m + enddo ! l + + ierr = NF_CLOSE(out_file_id) + + end Index: /BOL/Class_Reg/geo2reg.script =================================================================== --- /BOL/Class_Reg/geo2reg.script (revision 1585) +++ /BOL/Class_Reg/geo2reg.script (revision 1585) @@ -0,0 +1,35 @@ +#!/bin/bash + + +# +# localisation de l'executable: +repert='/u/fairhead/LMDZ/BOL/Class_Reg' +# +# liste des fichiers à traiter +liste_fichier='2L33_1m_2000*_histmth.nc' +# +# liste des variables a traiter +liste_variable='rneb rhum precip' + + +# + +cp $repert/config.def . + +var_list='aire,pourc_ter,w500' +for var in $liste_variable +do + var_list="$var_list,$var" +done + +ncrcat -c -v $var_list,t_ave_02592000 $liste_fichier toto.nc + +rm result.nc +for var in $liste_variable +do + $repert/geo2reg toto.nc $var + ncrcat -A out.nc result.nc > /dev/null 2>&1 +done +# +# un peu de menage +rm toto.nc Index: /BOL/IPCC_AR4/Makefile.linux =================================================================== --- /BOL/IPCC_AR4/Makefile.linux (revision 1585) +++ /BOL/IPCC_AR4/Makefile.linux (revision 1585) @@ -0,0 +1,16 @@ +F90=pgf90 +LIB_CMOR=/u/fairhead/IPCC/CMOR +#LIBNETCDF=-L/usr/local/pub/netcdf-3.5/lib -lnetcdf +LIBNETCDF=/u/fairhead/IPCC/CMOR/netcdf-3.5.1/lib +INCNETCDF=/u/fairhead/IPCC/CMOR/netcdf-3.5.1/include +LIBUDUNITS=-L/u/fairhead/IPCC/CMOR/udunits-1.11.7/src/lib/ -ludunits +LIBCDMS=-L/distrib/local/cdat-4.0b2_linux/lib/ -lcdms + +ts2IPCC: ts2IPCC.F90 + $(F90) -o ts2IPCC ts2IPCC.F90 -I$(LIB_CMOR) -L$(LIB_CMOR) -lcmor -I$(INCNETCDF) -L$(LIBNETCDF) -lnetcdf $(LIBUDUNITS) $(LIBCDMS) + +tostdlev: tostdlev.F90 plevel.o + $(F90) -o tostdlev tostdlev.F90 -I$(INCNETCDF) -L$(LIBNETCDF) -lnetcdf plevel.o + +plevel.o: plevel.F + $(F90) -c plevel.F Index: /BOL/IPCC_AR4/Makefile.sgi =================================================================== --- /BOL/IPCC_AR4/Makefile.sgi (revision 1585) +++ /BOL/IPCC_AR4/Makefile.sgi (revision 1585) @@ -0,0 +1,17 @@ +F90=f90 -n32 +LIB_CMOR=/home/rech/ces/rces011/CMOR +LIBIOIPSL=/home/rech/ces/rces011/ioipsl +#LIBNETCDF=-L/usr/local/pub/netcdf-3.5/lib -lnetcdf +LIBNETCDF=/home/rech/ces/rces011/netcdf-3.5.1/lib +INCNETCDF=/home/rech/ces/rces011/netcdf-3.5.1/include +LIBUDUNITS=-L/home/rech/ces/rces011/udunits-1.11.7/src/lib/ -ludunits +LIBCDMS=-L/home_b/rech/ces/rces599/PYTHON/CDAT/lib -lcdms + +ts2IPCC: ts2IPCC.F90 + $(F90) -o ts2IPCC ts2IPCC.F90 -I$(LIB_CMOR) -L$(LIB_CMOR) -lcmor -I$(INCNETCDF) -L$(LIBNETCDF) -lnetcdf $(LIBUDUNITS) $(LIBCDMS) + tostdlev: tostdlev.F90 plevel.o + $(F90) -o tostdlev tostdlev.F90 -I$(INCNETCDF) -L$(LIBNETCDF) -lnetcdf plevel.o + +plevel.o: plevel.F + $(F90) -c plevel.F + Index: /BOL/IPCC_AR4/config.def =================================================================== --- /BOL/IPCC_AR4/config.def (revision 1585) +++ /BOL/IPCC_AR4/config.def (revision 1585) @@ -0,0 +1,35 @@ +# Fichier de configuration pour le filtre IPSL2IPCC +# +# Parametres utilisés par cmor_setup: +# inpath = repertoire ou se trouve les tableaux MIP +inpath=../CMOR/Tables/IPCC/ +# file_action = que faire si un fichier à ecrire existe déjà +file_action = append +# verbosity = niveau de message +verbosity = 1 +# exit_control = que faire en cas d'erreur +exit_control = 2 +# +# Parametres utilisés par cmor_dataset +# repertoire = repertoire ou on range les fichier créés +repertoire=./Test +# experiment_ID = designation officielle de l'experience +experiment_ID = SRES A2 experiment +# institut = institut d'origine des données +institut = IPSL (Institut Pierre Simon Laplace, Paris, France) +# source = version du modele +source = IPSL-CM4_v1 +# realisation = numero de l'experience en cas d'ensemble +realisation = 1 +# contact = coordonnees de la personne responsable des donnees +contact = Sebastien Denvil, sebastien.denvil@ipsl.jussieu.fr +# hist_gen = historique general des donnees +hist_gen = YYYY/MM/JJ: data generated; YYYY/MM/JJ+1 data transformed +# comment = commentaires divers +comment = Test drive +# refs = references bibliographiques +refs = Dufresne et al, Journal of Climate, 2015, vol XX, p 136 +# +# Parametres utilises par cmor_variable +# hist_var = historique de la variable +hist_var = YYYY/MM/DD unit change Index: /BOL/IPCC_AR4/plevel.F =================================================================== --- /BOL/IPCC_AR4/plevel.F (revision 1585) +++ /BOL/IPCC_AR4/plevel.F (revision 1585) @@ -0,0 +1,126 @@ +c================================================================ +c================================================================ + SUBROUTINE plevel(ilon,ilev,lnew,pgcm,pres,Qgcm,Qpres) +c================================================================ +c================================================================ + + IMPLICIT none + + +c================================================================ +c +c Interpoler des champs 3-D u, v et g du modele a un niveau de +c pression donnee (pres) +c +c INPUT: ilon ----- nombre de points +c ilev ----- nombre de couches +c lnew ----- true si on doit reinitialiser les poids +c pgcm ----- pressions modeles +c pres ----- pression vers laquelle on interpolle +c Qgcm ----- champ GCM +c Qpres ---- champ interpolle au niveau pres +c +c================================================================ +c +c arguments : +c ----------- + + INTEGER ilon, ilev + logical lnew + + REAL pgcm(ilon,ilev) + REAL Qgcm(ilon,ilev) + real pres + REAL Qpres(ilon) + +c local : +c ------- + +c INTEGER lt(ilon), lb(ilon) +c REAL ptop, pbot, aist(ilon), aisb(ilon) + REAL ptop, pbot + + integer, dimension(:), allocatable :: lt, lb + real, dimension(:), allocatable :: aist, aisb + save lt,lb,ptop,pbot,aist,aisb + + INTEGER i, k + logical first + SAVE first + data first /.true./ +c + + if (first) THEN + allocate (lt(ilon)) + allocate (lb(ilon)) + allocate (aist(ilon)) + allocate (aisb(ilon)) + first=.false. + endif +c===================================================================== + if (lnew) THEN + +c on réinitialise les réindicages et les poids +c===================================================================== + + +c Chercher les 2 couches les plus proches du niveau a obtenir +c +c Eventuellement, faire l'extrapolation a partir des deux couches +c les plus basses ou les deux couches les plus hautes: + DO 130 i = 1, ilon + IF ( ABS(pres-pgcm(i,ilev) ) .LT. + . ABS(pres-pgcm(i,1)) ) THEN + lt(i) = ilev ! 2 + lb(i) = ilev-1 ! 1 + ELSE + lt(i) = 2 + lb(i) = 1 + ENDIF + 130 CONTINUE + DO 150 k = 1, ilev-1 + DO 140 i = 1, ilon + pbot = pgcm(i,k) + ptop = pgcm(i,k+1) + IF (ptop.LE.pres .AND. pbot.GE.pres) THEN + lt(i) = k+1 + lb(i) = k + ENDIF + 140 CONTINUE + 150 CONTINUE +c +c Interpolation lineaire: +c + DO i = 1, ilon +c interpolation en logarithme de pression: +c +c ... Modif . P. Le Van ( 20/01/98) .... +c Modif Frédéric Hourdin (3/01/02) + + aist(i) = LOG( pgcm(i,lb(i))/ pres ) + . / LOG( pgcm(i,lb(i))/ pgcm(i,lt(i)) ) + aisb(i) = LOG( pres / pgcm(i,lt(i)) ) + . / LOG( pgcm(i,lb(i))/ pgcm(i,lt(i))) + enddo + + endif ! lnew + +c====================================================================== +c inteprollation +c====================================================================== + + do i=1,ilon + Qpres(i)= Qgcm(i,lb(i))*aisb(i)+Qgcm(i,lt(i))*aist(i) + enddo +c +c Je mets les vents a zero quand je rencontre une montagne + do i = 1, ilon + if (pgcm(i,1).LT.pres) THEN +c Qpres(i)=1e33 + Qpres(i)=1e+20 + endif + enddo + +c + RETURN + END Index: /BOL/IPCC_AR4/table.def =================================================================== --- /BOL/IPCC_AR4/table.def (revision 1585) +++ /BOL/IPCC_AR4/table.def (revision 1585) @@ -0,0 +1,70 @@ +# +# tableau de correspondance entre les noms des variables du modele +# ipsl et des variables IPCC +# +# Nom IPSL | unites IPSL | sens | nom IPCC | table IPCC +# |positif| +slp | Pa | -- | psl | IPCC_table_A1 +t2m | K | -- | tas | IPCC_table_A1 +precip | kg m-2 s-1 | -- | pr | IPCC_table_A1 + | | -- | mrsos | IPCC_table_A1 + | | -- | mrso | IPCC_table_A1 + | | -- | tauu | IPCC_table_A1 + | | -- | tauv | IPCC_table_A1 +snow_mass | kg m-2 | -- | snd | IPCC_table_A1 + | | -- | hfls | IPCC_table_A1 +sens | W m-2 | -- | hfss | IPCC_table_A1 +LWdnSFC | W m-2 | -- | rlds | IPCC_table_A1 +LWupSFC | W m-2 | -- | rlus | IPCC_table_A1 +SWdnSFC | W m-2 | -- | rsds | IPCC_table_A1 +SWupSFC | W m-2 | up | rsus | IPCC_table_A1 +tsol | K | -- | ts | IPCC_table_A1 +psol | Pa | -- | ps | IPCC_table_A1 +snow | kg m-2 s-1 | -- | prsn | IPCC_table_A1 +pluc | kg m-2 s-1 | -- | prc | IPCC_table_A1 + | | -- | prw | IPCC_table_A1 + | | -- | mrfso | IPCC_table_A1 + | | -- | mrros | IPCC_table_A1 + | | -- | mrro | IPCC_table_A1 + | | -- | snw | IPCC_table_A1 + | | -- | snc | IPCC_table_A1 + | | -- | snm | IPCC_table_A1 +u10m | m s-1 | -- | uas | IPCC_table_A1 +v10m | m s-1 | -- | vas | IPCC_table_A1 +q2m | kg kg-1 | -- | huss | IPCC_table_A1 +SWdnTOA | W m-2 | -- | rsdt | IPCC_table_A1 +SWupTOA | W m-2 | -- | rsut | IPCC_table_A1 +topl | W m-2 | -- | rlut | IPCC_table_A1 + | | -- | rtmt | IPCC_table_A1 +SWdn200 | W m-2 | -- | rsntp | IPCC_table_A1 +LWup200 | W m-2 | -- | rlntp | IPCC_table_A1 +SWdn200clr | W m-2 | -- | rsntpcs | IPCC_table_A1 +LWup200clr | W m-2 | -- | rlntpcs | IPCC_table_A1 +SWdnSFCclr | W m-2 | -- | rsdscs | IPCC_table_A1 +SWupSFCclr | W m-2 | -- | rsuscs | IPCC_table_A1 +LWdnSFCclr | W m-2 | -- | rldscs | IPCC_table_A1 + | | -- | rlutcs | IPCC_table_A1 +SWupTOAclr | W m-2 | -- | rsutcs | IPCC_table_A1 +cldt | | -- | clt | IPCC_table_A1 +lwp | kg m-2 | -- | clwvi | IPCC_table_A1 +iwp | kg m-2 | -- | clivi | IPCC_table_A1 + | | -- | orog | IPCC_table_A1 +pourc_ter | | -- | sftlf | IPCC_table_A1 +pourc_lic | | -- | sftgif | IPCC_table_A1 + | | -- | mrsofc | IPCC_table_A1 +rneb | % | -- | cl | IPCC_table_A1 +temp | K | -- | ta | IPCC_table_A1 +vitu | m s-1 | -- | ua | IPCC_table_A1 +vitv | m s-1 | -- | va | IPCC_table_A1 +ovap | kg kg-1 | -- | hus | IPCC_table_A1 + | | -- | wap | IPCC_table_A1 +geop | m | -- | zg | IPCC_table_A1 +rhum | | -- | hur | IPCC_table_A1 +ozone | | -- | tro3 | IPCC_table_A1 + | | -- | clisccp | IPCC_table_A1 + | | -- | trsul | IPCC_table_A1 + | | -- | trsult | IPCC_table_A1 + | | -- | prveg | IPCC_table_A1 + | | -- | evspsblveg | IPCC_table_A1 + | | -- | zmla | IPCC_table_A1 + | | -- | rootd | IPCC_table_A1 Index: /BOL/IPCC_AR4/tostdlev.F90 =================================================================== --- /BOL/IPCC_AR4/tostdlev.F90 (revision 1585) +++ /BOL/IPCC_AR4/tostdlev.F90 (revision 1585) @@ -0,0 +1,247 @@ + + PROGRAM tostdlev + +! +! Pour passer a posteriori des champs 3d du modele sur des niveaux standards +! +! L. Fairhead 2004/12 +! +! Ce programme est appelé avec un argument, le nom du fichier à traiter +! Il nécessite aussi un fichier config.def contenant diverses informations +! (voir plus bas) et l'accès à un tableau faisant la correspondance entre +! les noms de variables du modèle et les noms imposés par l'IPCC +! Pour l'instant on ne traite que les fichiers contenant la serie +! temporelle d'une seule variable + + implicit none + +#include "netcdf.inc" + + INTEGER :: lunout, ierr, icount, ivar, nvars, nlev, len + integer :: press_id, data_id, newdata_id + INTEGER :: lonid, ilon, newlonid + INTEGER :: latid, ilat, newlatid + INTEGER :: levid, ilev, newlevid + INTEGER :: vartype, ndims + INTEGER :: timid, itime, newtimid + INTEGER :: varlonid, newvarlonid + INTEGER :: varlatid, newvarlatid + INTEGER :: vartimid, newvartimid + INTEGER :: varlevid, newvarlevid + INTEGER :: varid, newvarid + integer :: presvar_id + INTEGER, dimension(4) :: dimids + integer :: natts + CHARACTER (len=80) :: varname, attname + INTEGER, DIMENSION(4) :: start, count + + REAL, DIMENSION(:), ALLOCATABLE :: lon, lat, time + REAL, DIMENSION(:,:,:), ALLOCATABLE :: pression, champ, newchamp + + +! niveaux standards: + INTEGER,parameter :: nlevstd = 17 + real :: rlevstd(nlevstd) + DATA rlevstd /100000., 92500., 85000., 70000., & + & 60000., 50000., 40000., 30000., 25000., 20000., & + & 15000., 10000., 7000., 5000., 3000., 2000., 1000./ +! +! quelques initialisations + lunout = 6 + +! +! Ouverture du fichier contenant les pressions modele + ierr = nf_open('pression.nc', NF_NOWRITE, press_id) + IF (ierr /= NF_NOERR) then + WRITE(lunout,*)NF_STRERROR(ierr) + stop + endif + +! +! Ouverture du fichier contenant le champ a interpoler + ierr = nf_open('data.nc', NF_NOWRITE, data_id) + IF (ierr /= NF_NOERR) then + WRITE(lunout,*)NF_STRERROR(ierr) + stop + ENDIF + +! +! Ouverture du fichier contenant le champ interpolé + ierr = nf_create('newdata.nc', NF_CLOBBER, newdata_id) + IF (ierr /= NF_NOERR) then + WRITE(lunout,*)NF_STRERROR(ierr) + stop + endif + +! +! Début lecture fichier origine/ecriture fichier sortie +! +! Definition des dimensions: + ierr = nf_inq_ndims(data_id,ndims) + DO icount = 1, ndims + ierr = nf_inq_dim(data_id, icount, varname, len) + select case (trim(varname)) + case ('lon') + ilon = len +! Definition de la longitude + ierr = nf_def_dim(newdata_id,varname,ilon, newlonid) + case ('lat') + ilat = len +! Definition de la latitude + ierr = nf_def_dim(newdata_id,varname,ilat, newlatid) + case('presnivs') + ilev = len +! Definition niveaux verticaux + ierr = nf_def_dim(newdata_id,'presnivs',nlevstd,newlevid) + case ('time_counter') + itime = len +! Definition du temps + ierr = nf_def_dim(newdata_id,varname,itime, newtimid) + case default + WRITE(lunout,*)'je ne reconnais pas cette dimension: ',varname + stop + end select + enddo +! +! Definition des variables +! +! +! + ierr = nf_inq_nvars(data_id, nvars) + DO ivar = 1, nvars + ierr = nf_inq_var(data_id, ivar, varname, vartype, ndims, dimids, natts) + if (ierr /= 0) call handle_err(ierr) + selectcase (trim(varname)) + case('lon') +! definition de la longitude + ierr = nf_inq_varid(data_id,varname, varlonid) + dimids(1)=newlonid + ierr = nf_def_var(newdata_id, varname, vartype , ndims, dimids,newvarlonid) +! recopie des attributs de la variable: + DO icount = 1, natts + ierr = nf_inq_attname(data_id, varlonid, icount, attname) + ierr = nf_copy_att(data_id, varlonid, attname, newdata_id, newvarlonid) + enddo + case('lat') +! definition de la latitude + ierr = nf_inq_varid(data_id, varname, varlatid) + dimids(1)=newlatid + ierr = nf_def_var(newdata_id, varname, vartype , ndims, dimids,newvarlatid) +! recopie des attributs de la variable: + DO icount = 1, natts + ierr = nf_inq_attname(data_id, varlatid, icount, attname) + ierr = nf_copy_att(data_id, varlatid, attname, newdata_id, newvarlatid) + enddo + case('time_counter') +! definition du temps + ierr = nf_inq_varid(data_id, varname, vartimid) + dimids(1)=newtimid + ierr = nf_def_var(newdata_id, varname, vartype, ndims, dimids,newvartimid) +! recopie des attributs de la variable: + DO icount = 1, natts + ierr = nf_inq_attname(data_id, vartimid, icount, attname) + ierr = nf_copy_att(data_id, vartimid, attname, newdata_id, newvartimid) + enddo + case('presnivs') +! definition des niveaux de pression + ierr = nf_inq_varid(data_id,varname,varlevid) + dimids(1)=newlevid + ierr = nf_def_var(newdata_id, varname, vartype, ndims, dimids,newvarlevid) +! recopie des attributs de la variable: + ierr = nf_copy_att(data_id, varlevid, 'units', newdata_id, newvarlevid) + ierr = nf_copy_att(data_id, varlevid, 'title', newdata_id, newvarlevid) + ierr = nf_copy_att(data_id, varlevid, 'long_name', newdata_id, newvarlevid) + case default +! normalement il ne reste que la variable à interpoler 4d + IF (ndims /= 4) then + WRITE(lunout,*)'La variable principale du fichier n''est pas 4D' + stop + endif + ierr = nf_inq_varid(data_id,varname,varid) + if (ierr /= 0) call handle_err(ierr) + dimids(1) = newlonid + dimids(2) = newlatid + dimids(3) = newlevid + dimids(4) = newtimid + ierr = nf_def_var(newdata_id, varname, vartype, ndims, dimids,newvarid) + if (ierr /= 0) call handle_err(ierr) + DO icount = 1, natts + ierr = nf_inq_attname(data_id, varid, icount, attname) + ierr = nf_copy_att(data_id, varid, attname, newdata_id, newvarid) + enddo + end select + enddo +! + +! +! fermeture du mode definition du fichier + ierr = nf_enddef(newdata_id) +! +! ecriture des variables: +! longitude + allocate(lon(ilon)) + ierr = nf_get_var_real(data_id, varlonid, lon) + ierr = nf_put_var_real(newdata_id, newvarlonid,lon) +! +! latitude + allocate(lat(ilat)) + ierr = nf_get_var_real(data_id, varlatid, lat) + ierr = nf_put_var_real(newdata_id, newvarlatid, lat) +! +! niveaux de pression + ierr = nf_put_var_real(newdata_id, newvarlevid, rlevstd) +! +! temps + allocate(time(itime)) + ierr = nf_get_var_real(data_id, vartimid, time) + ierr = nf_put_var_real(newdata_id, newvartimid, time) + +! Traitement de la variable + ALLOCATE(pression(ilon,ilat,ilev)) + ierr = nf_inq_varid(press_id, 'pres', presvar_id) + ALLOCATE(champ(ilon,ilat,ilev)) + ALLOCATE(newchamp(ilon,ilat,nlevstd)) + start(1) = 1 + start(2) = 1 + start(3) = 1 + COUNT(1) = ilon + COUNT(2) = ilat + COUNT(3) = ilev + count(4) = 1 + DO icount = 1, itime + COUNT(3) = ilev + start(4) = icount + ierr = nf_get_vara_real(press_id, presvar_id, start, count, pression) + if (ierr /= 0) call handle_err(ierr) + ierr = nf_get_vara_real(data_id, varid, start, count, champ) + if (ierr /= 0) call handle_err(ierr) + DO nlev = 1, nlevstd + CALL plevel(ilon*ilat, ilev,.TRUE., pression, rlevstd(nlev),champ, & + & newchamp(:,:,nlev)) + enddo + count(3) = nlevstd + ierr = nf_put_vara_real(newdata_id, newvarid, start, count, newchamp) + if (ierr /= 0) call handle_err(ierr) + enddo + +! On ferme! + ierr = nf_close(press_id) + ierr = nf_close(data_id) + ierr = nf_close(newdata_id) + +contains + +subroutine handle_err(status) + use netcdf + + implicit none + integer, intent(in) :: status + + write(lunout,*)nf90_strerror(status) + stop + +end subroutine handle_err + + +END PROGRAM tostdlev + Index: /BOL/IPCC_AR4/ts2IPCC.F90 =================================================================== --- /BOL/IPCC_AR4/ts2IPCC.F90 (revision 1585) +++ /BOL/IPCC_AR4/ts2IPCC.F90 (revision 1585) @@ -0,0 +1,433 @@ +! +! $Header$ +! +PROGRAM ts2IPCC + +! +! Filtre permettant de transformer les fichiers de serie temporelle a une +! variable de l'IPSL en fichiers acceptables par le PCMDI/IPCC. +! +! Utilisation de la bibliothèque CMOR du PCMDI +! +! L. Fairhead 2004/08 +! +! Ce programme est appelé avec un argument, le nom du fichier à traiter +! Il nécessite aussi un fichier config.def contenant diverses informations +! (voir plus bas) et l'accès à un tableau faisant la correspondance entre +! les noms de variables du modèle et les noms imposés par l'IPCC +! Pour l'instant on ne traite que les fichiers contenant la serie +! temporelle d'une seule variable + + use cmor_users_functions + use netcdf + + implicit none + +#include "netcdf.inc" + + CHARACTER (len=256) :: orig_file ! nom du fichier à traiter + character (len=512) :: line_read + CHARACTER (len=128) :: inpath, contact, hist_gen,repert,instit + CHARACTER (len=128) :: hist_var,expt_id,source,comment,refs + CHARACTER (len=20) :: action + character (len=20) :: first_part + character (len=1004) :: second_part + CHARACTER (len=20), DIMENSION(100) :: ipsl_name, ipsl_units, ipsl_pos + CHARACTER (len=20), DIMENSION(100) :: ipcc_name, ipcc_table + CHARACTER (len=80) :: varname, units, namedim + + INTEGER :: orig_file_id, nvars, ndims + INTEGER :: verbos, exit_ctl, realis, indice,index_table + INTEGER :: iargc, iostat, ierr + INTEGER :: i,idim + INTEGER, ALLOCATABLE, DIMENSION(:) :: dimids,axis_ids,lendim + INTEGER :: latid, lonid, vertid, timeid + INTEGER :: varid, cmorvarid + INTEGER :: ilat, ilon, ivert, itime + INTEGER :: lunout ! device de sortie + + logical :: found = .false. + + REAL, ALLOCATABLE, DIMENSION(:) :: lon, lat, vert, time + REAL, ALLOCATABLE, DIMENSION(:) :: lon_bounds, lat_bounds + REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: donnees + DOUBLE PRECISION, DIMENSION(1) :: rdate + real :: missing_value + + external iargc + +! +! quelques initialisations + lunout = 6 + varname = 'xxxxxxxx' + +! +! On vérifie que l'appel au programme a bien un argument: + CALL getarg(1, orig_file) + IF (iargc() == 0 .OR. orig_file == '-h') then + WRITE(lunout,*)' ' + WRITE(lunout,*)' Utilisation de ce programme: ' + WRITE(lunout,*)' ./ts2IPCC nom_de_fichier [variable]' + WRITE(lunout,*)' ou nom_de_fichier est le nom du fichier a traiter' + WRITE(lunout,*)' et variable la variable a traiter [optionel]' + WRITE(lunout,*)' ' + WRITE(lunout,*)' ./ts2IPCC -h sort ce message' + WRITE(lunout,*)' ' + stop + ENDIF + if (iargc() == 2) then + CALL getarg(2, varname) + endif + +! +! Lecture du fichier de configuration + OPEN (20, IOSTAT=iostat, file='config.def',form='formatted') + IF (iostat /= 0) then + WRITE(lunout,*)'Erreur ouverture du fichier config.def' + stop + endif + + config: do + READ(20,'(A)',iostat=iostat)line_read + IF (iostat /= 0) exit + line_read = TRIM(line_read) + IF (INDEX(line_read, '#') /= 1) THEN + first_part = trim(line_read(1:INDEX(line_read, '=')-1)) + second_part = trim(line_read(INDEX(line_read, '=')+1:)) + selectcase(first_part) + case('inpath') + inpath = trim(second_part) + case('file_action') + action = trim(second_part) + case('verbosity') + READ(second_part,'(i)') verbos + case('exit_control') + READ(second_part,'(i)') exit_ctl + case('repertoire') + repert = trim(second_part) + case('experiment_ID') + expt_id = trim(second_part) + case('institut') + instit = trim(second_part) + case('source') + source = trim(second_part) + case('realisation') + READ(second_part,'(i)') realis + case('hist_gen') + hist_gen = trim(second_part) + case('comment') + comment = trim(second_part) + case('refs') + refs = trim(second_part) + case('hist_var') + hist_var = trim(second_part) + case('contact') + contact = trim(second_part) + end select + endif + enddo config + if (iostat > 0) then + WRITE(lunout,*)'Probleme de lecture du fichier config.def, iostat = ',iostat + stop + endif + close(20) + +! +! Lecture du tableau de correspondance nom IPSL <=> nom IPCC + OPEN (20, IOSTAT=iostat, file='table.def',form='formatted') + IF (iostat /= 0) then + WRITE(lunout,*)'Erreur ouverture du fichier table.def' + stop + endif + indice = 0 + table: do + READ(20,'(A)',iostat=iostat)line_read + IF (iostat /= 0) exit + line_read = TRIM(line_read) + IF (INDEX(line_read, '#') /= 1) THEN + indice = indice + 1 + ipsl_name(indice) = trim(line_read(1:INDEX(line_read, '|')-1)) + line_read = trim(line_read(INDEX(line_read, '|')+1:)) + ipsl_units(indice) = trim(line_read(1:INDEX(line_read, '|')-1)) + line_read = trim(line_read(INDEX(line_read, '|')+1:)) + ipsl_pos(indice) = trim(line_read(1:INDEX(line_read, '|')-1)) + line_read = trim(line_read(INDEX(line_read, '|')+1:)) + ipcc_name(indice) = trim(line_read(1:INDEX(line_read, '|')-1)) + ipcc_table(indice) = trim(line_read(INDEX(line_read, '|')+1:)) + endif + enddo table + indice = indice - 1 + close(20) +! DO i = 1, indice +! WRITE(lunout,*)ipsl_name(i),ipsl_units(i),ipcc_name(i),ipcc_table(i) +! enddo + +! +! Ouverture du fichier a traiter + ierr = nf_open(orig_file, NF_NOWRITE, orig_file_id) + IF (ierr /= NF_NOERR) then + WRITE(lunout,*)NF_STRERROR(ierr) + stop + endif +! +! trouver la variable a traiter, c'est une variable a 3 ou 4 dimensions + ierr = nf_inq_nvars(orig_file_id, nvars) + IF (ierr /= NF_NOERR) then + WRITE(lunout,*)NF_STRERROR(ierr) + stop + endif + + i = 0 + if (varname == 'xxxxxxxx') then + DO while (.not.found) + i = i + 1 + if (i > nvars) then + WRITE(lunout,*)' pas de variable 3d ou 4d trouvee' + stop + endif + ierr = nf_inq_varname(orig_file_id, i, varname) + IF (ierr /= NF_NOERR) then + WRITE(lunout,*)NF_STRERROR(ierr) + stop + endif + ierr = nf_inq_varndims(orig_file_id, i, ndims) + IF (ierr /= NF_NOERR) then + WRITE(lunout,*)NF_STRERROR(ierr) + stop + endif + if (ndims > 2) found = .true. + enddo + else + ierr = nf_inq_varid(orig_file_id, varname, varid) + IF (ierr /= NF_NOERR) then + WRITE(lunout,*)NF_STRERROR(ierr) + stop + endif + endif + +! +! recherche de la correspondance nom IPSL <=> nom IPCC + found = .false. + i = 0 + do while (.not. found) + i = i + 1 + if (i > indice) then + WRITE(lunout,*)'La variable ',trim(varname),' n''est pas dans le tableau de correspondance table.def' + stop + endif + IF (varname == ipsl_name(i)) THEN + index_table = i + found = .true. + endif + enddo + + WRITE(lunout,*)' found variable = ', trim(varname) + WRITE(lunout,*)' ipcc_name = ', trim(ipcc_name(index_table)) + + +! +! Initialisation CMOR + ierr = cmor_setup(inpath=inpath, netcdf_file_action=action,set_verbosity=verbos,& + & exit_control=exit_ctl) + IF (ierr /= 0) then + WRITE(lunout,*)'Probleme dans cmor_setup, ierr = ', ierr + endif + +! +! Initialisation dataset + ierr = cmor_dataset(outpath=repert, & + & experiment_id=expt_id, & + & institution=instit, & + & source=source, & + & calendar='360_day', & + & realization=realis, & + & contact=contact, & + & history=hist_gen, & + & comment=comment, & + & references=refs) + IF (ierr /= 0) then + WRITE(lunout,*)'Probleme dans cmor_dataset, ierr = ', ierr + endif + +! +! Definition des axes + + ierr = nf90_inq_varid(orig_file_id,TRIM(varname), varid) + if (ierr /= 0) call handle_err(ierr) + ierr = nf90_Inquire_Variable(orig_file_id, varid, ndims = ndims) + if (ierr /= 0) call handle_err(ierr) + allocate (dimids(ndims)) + allocate (axis_ids(ndims)) + allocate (lendim(ndims)) + ierr = nf90_Inquire_Variable(orig_file_id, varid, dimids = dimids) + if (ierr /= 0) call handle_err(ierr) + + do idim = 1, ndims + ierr = nf90_Inquire_Dimension(orig_file_id, dimids(idim), & + name = namedim, len = lendim(idim)) + if (ierr /= 0) call handle_err(ierr) + units=' ' + selectcase(trim(namedim)) + case('lat') +! lecture de la latitude: + allocate(lat(lendim(idim))) + ierr = nf_inq_varid(orig_file_id, namedim, latid) + if (ierr /= 0) call handle_err(ierr) + ierr = nf_get_var_real(orig_file_id, latid, lat) + if (ierr /= 0) call handle_err(ierr) + ierr = nf_get_att_text(orig_file_id, latid, 'units', units) + if (ierr /= 0) call handle_err(ierr) + allocate(lat_bounds(lendim(idim)+1)) + DO i = 2, lendim(idim) + lat_bounds(i) = lat(i-1) - (lat(i-1) - lat(i))/2 + enddo + lat_bounds(1) = lat(1) + lat_bounds(lendim(idim)+1) = lat(lendim(idim)) +! definition de la latitude + axis_ids(idim) = cmor_axis( & + table=trim(ipcc_table(index_table)), & + table_entry='latitude', & + units=trim(units), & + length=lendim(idim), & + coord_vals=lat, & + cell_bounds=lat_bounds) +! +! + case('lon') +! lecture de la longitude: + allocate(lon(lendim(idim))) + ierr = nf_inq_varid(orig_file_id, namedim, lonid) + if (ierr /= 0) call handle_err(ierr) + ierr = nf_get_var_real(orig_file_id, lonid, lon) + if (ierr /= 0) call handle_err(ierr) + ierr = nf_get_att_text(orig_file_id, lonid, 'units', units) + if (ierr /= 0) call handle_err(ierr) + ALLOCATE(lon_bounds(lendim(idim)+1)) + DO i = 2, lendim(idim) + lon_bounds(i) = lon(i-1) - (lon(i-1) - lon(i))/2 + enddo + lon_bounds(1) = lon(1) - (lon_bounds(3) -lon_bounds(2))/2. + lon_bounds(lendim(idim)+1) = lon(lendim(idim)) + (lon_bounds(lendim(idim))-lon_bounds(lendim(idim)-1))/2. + +! definition de la longitude + axis_ids(idim) = cmor_axis( & + table=trim(ipcc_table(index_table)), & + table_entry='longitude', & + units=trim(units), & + length=lendim(idim), & + coord_vals=lon, & + cell_bounds=lon_bounds) +! +! + case('presnivs') +! lecture de la verticale: + allocate(vert(lendim(idim))) + ierr = nf_inq_varid(orig_file_id, namedim, vertid) + if (ierr /= 0) call handle_err(ierr) + ierr = nf_get_var_real(orig_file_id, vertid, vert) + if (ierr /= 0) call handle_err(ierr) + ierr = nf_get_att_text(orig_file_id, vertid, 'units', units) + if (ierr /= 0) call handle_err(ierr) +! +! definition de la verticale + if (units == 'mb') units='Pa' + axis_ids(idim) = cmor_axis( & + table=trim(ipcc_table(index_table)), & + table_entry='pressure', & + units=trim(units), & + length=lendim(idim), & + coord_vals=vert) +! +! + case('time_counter') +! definition du temps + if (idim /= ndims) then + write(lunout,*)'la dimension temps doit etre la derniere dimension' + stop + endif + allocate(time(lendim(idim))) + ierr = nf_inq_varid(orig_file_id,namedim,timeid) + if (ierr /=0) call handle_err(ierr) + ierr = nf_get_var_real(orig_file_id, timeid, time) + if (ierr /=0) call handle_err(ierr) + ierr = nf_get_att_text(orig_file_id,timeid, 'units', units) + if (ierr /=0) call handle_err(ierr) + axis_ids(idim) = cmor_axis( & + table=trim(ipcc_table(index_table)), & + table_entry='time', & + units=trim(units), & + length=lendim(idim), & + interval='30 minutes') + itime= lendim(idim) + case default + write(lunout,*)'Dimension: ', trim(namedim),' non reconnue' + stop + endselect + enddo + +! +! Definition de la variable a ecrire + units=' ' + ierr = nf_inq_varid(orig_file_id,TRIM(varname), varid) + if (ierr /= 0) call handle_err(ierr) + ierr = nf_get_att_text(orig_file_id, varid, 'units', units) + if (ierr /= 0) call handle_err(ierr) + ierr = nf_get_att_real(orig_file_id, varid, 'missing_value', missing_value) + if (ierr /= 0) call handle_err(ierr) + cmorvarid = cmor_variable( & + table=trim(ipcc_table(index_table)), & + table_entry=trim(ipcc_name(index_table)), & + units=trim(units), & + axis_ids= axis_ids, & + missing_value=real(missing_value), & + positive = trim(ipsl_pos(index_table)), & + original_name=trim(varname)) +! +! Lecture de la variable + if (ndims == 3) then + ALLOCATE (donnees(lendim(1), lendim(2), 1, lendim(3) )) + else if (ndims ==4) then + allocate (donnees(lendim(1), lendim(2), lendim(3), lendim(4) )) + endif + ierr = nf_get_var_real(orig_file_id, varid, donnees) +! +! Ecriture de la variable + + DO i = 1, itime + rdate(1) = dble(time(i)) + ierr = cmor_write( & + var_id = cmorvarid, & + DATA = REAL(donnees(:,:,:,i)), & + ntimes_passed = 1, & + time_vals = rdate) + enddo +! +! Fin CMOR + ierr = cmor_close() + IF (ierr /= 0) then + WRITE(lunout,*)'Probleme dans cmor_close, ierr = ', ierr + endif + +! +! fermeture fichier originel + ierr = nf_close(orig_file_id) + IF (ierr /= NF_NOERR) then + WRITE(lunout,*)NF_STRERROR(ierr) + stop + endif + +contains + +subroutine handle_err(status) + use netcdf + + implicit none + integer, intent(in) :: status + + write(lunout,*)nf90_strerror(status) + stop + +end subroutine handle_err + +END PROGRAM ts2IPCC +