!> \file out_cptr_mod.f90 !! Module avec les routines d'ecriture et de lecture de fichiers cptr !< !> \namespace out_cptr !! This module gathers routines to read and write cptr files !! \author Vincent Peyaud !! \date janvier 2006 !! @note attention : en ecriture fichier runname.cptr !! @note en lecture fichier runname.CPTR !! @note il faut renommer le fichier ! !< !---------------------------------------------------------- !> subroutine: read_recovery !!Reprise d'un fichier compteur !! Used module: - use module3D_phy !! - use netcdf !! - use io_netcdf_grisli !! - use tracer_vars !! ! !> subroutine read_recovery(icpt) !---- reprise d'un fichier compteur ----- use module3D_phy use netcdf use io_netcdf_grisli use tracer_vars !aurel neem implicit none character(len=20) :: titre integer :: nzzm integer :: icpt ! icompteur 1-> tout (topo, T, U) ! 2 -> tout sauf topo ! 3 -> tout sauf topo et eau basale real,dimension(nx,ny) :: bidon !< sert à sauter les lectures de Ss,H,B,... ! Variables specifiques pour les fichier .nc real*8, dimension(:,:), pointer :: tab => null() !< tableau 2d real ecrit dans le fichier real*8, dimension(:,:,:), pointer :: tab1 => null() !< tableau 3d real real*8, dimension(:,:,:), pointer :: tab1T => null() !< tableau 3d real pour la temperature real*8, dimension(:), pointer :: liste_time => null() ! rappel de f90 ! reprcptr est un string (ici le nom du fichier recovery) ! index (reprcptr,'.cptr') renvoie : ! 0, si 'cptr' n'est pas dans reprcptr ! la position de la premiere occurence de cptr sinon if (ncdf_type.eq.0) call lect_netcdf_type !< pour lire la valeur de netcdf_type (machine dependant) ascii_nc: if ((index(reprcptr,'.cptr').ne.0) .or.(index(reprcptr,'.CPTR').ne.0)) then !reprise fichier cptr if (itracebug.eq.1) call tracebug(' Entree dans routine read_recovery : ascii') open(num_forc,file=trim(reprcptr)) read(num_forc,*) time ! temps de la reprise !----------------------- if (tbegin.gt.1.d9) then ! si tbegin est tres grand , on reprend le temps du cptr tbegin=time else time=tbegin endif read(num_forc,*) titre read(num_forc,*) ni,nzz,nzzm,nxx,nyy read(num_forc,*) if ((nxx.ne.nx).or.(nyy.ne.ny).or.(nzz.ne.nz) & .or.(nzzm.ne.nzm)) write(6,*) 'attention', & 'les tailles de tableaux ne sont pas compatibles' ! lecture des tableaux 2D read(num_forc,*) read(num_forc,*) lect_topo: if (icpt.eq.1) then ! lecture topo seulement si icompteur = 1 read(num_forc,*) S read(num_forc,*) read(num_forc,*) H read(num_forc,*) read(num_forc,*) B read(num_forc,*) read(num_forc,*) Bsoc read(num_forc,*) else read(num_forc,*) bidon read(num_forc,*) read(num_forc,*) bidon read(num_forc,*) read(num_forc,*) bidon read(num_forc,*) read(num_forc,*) bidon read(num_forc,*) end if lect_topo ! lecture temperature dans tous les cas read(num_forc,*) Ibase read(num_forc,*) read(num_forc,*) bmelt read(num_forc,*) lect_Hwater: if (icpt.ne.3) then ! pas de lecture Hwater si icompteur = 3 read(num_forc,*) Hwater read(num_forc,*) else read(num_forc,*) bidon read(num_forc,*) Hwater(:,:)=0. end if lect_Hwater ! lecture vitesses dans tous les cas. read(num_forc,*) Uxbar read(num_forc,*) read(num_forc,*) Uybar ! tableaux 3d read(num_forc,*) read(num_forc,*) read(num_forc,*,end=22) T 22 continue read(num_forc,*,end=23) read(num_forc,*,end=23) read(num_forc,*,end=23) Ux 23 continue read(num_forc,*,end=24) read(num_forc,*,end=24) read(num_forc,*,end=24) Uy 24 continue ! aurel pour faire des cptr avec les traceurs : read(num_forc,*,end=25) read(num_forc,*,end=25) read(num_forc,*,end=25) Xdep_out !xdepk ou xdep ? a voir 25 continue read(num_forc,*,end=26) read(num_forc,*,end=26) read(num_forc,*,end=26) Ydep_out 26 continue read(num_forc,*,end=27) read(num_forc,*,end=27) read(num_forc,*,end=27) tdep_out 27 continue tdep_out(:,:,:)=tdep_out(:,:,:)+time close(num_forc) else if ((index(reprcptr,'.nc').ne.0)) then ! Modif hassine pour la reprise d'un fichier .nc if (itracebug.eq.1) call tracebug(' Entree dans routine read_recovery : netcdf') if (.not.associated(liste_time)) then allocate(liste_time(1)) end if if (.not.associated(tab)) then allocate(tab(nx,ny),tab1(nx,ny,nz),tab1T(nx,ny,nz+nzm)) end if call Read_ncdf_var('temps',TRIM(reprcptr),liste_time) ! temps de la reprise !----------------------- if (tbegin.gt.1.d9) then ! si tbegin est tres grand , on reprend le temps du .nc tbegin=liste_time(1) time = tbegin else time = tbegin endif if (itracebug.eq.1) write(num_tracebug,*) tbegin, liste_time(1) ! TBEGIN=liste_time(1)-liste_time(1) call Read_ncdf_dim('x' ,TRIM(reprcptr),nxx) call Read_ncdf_dim('y' ,TRIM(reprcptr),nyy) call Read_ncdf_dim('z' ,TRIM(reprcptr),nzz) call Read_ncdf_dim('zm',TRIM(reprcptr),nzzm) if ((nxx.ne.nx).or.(nyy.ne.ny).or.(nzz.ne.nz) & .or.(nzzm-nzz.ne.nzm)) write(6,*) 'attention', & 'les tailles de tableaux ne sont pas compatibles' ! lecture des tableaux 2D ! ----------------------------- if (icpt.eq.1) then ! lecture topo seulement si icompteur = 1 call Read_ncdf_var('S',TRIM(reprcptr),tab) S(:,:)=tab(:,:) call Read_ncdf_var('H',TRIM(reprcptr),tab) H(:,:)=tab(:,:) call Read_ncdf_var('B',TRIM(reprcptr),tab) B(:,:)=tab(:,:) call Read_ncdf_var('BSOC',TRIM(reprcptr),tab) Bsoc(:,:)=tab(:,:) end if ! lecture temperature dans tous les cas call Read_ncdf_var('IBASE',TRIM(reprcptr),tab) Ibase(:,:)=tab(:,:) call Read_ncdf_var('BMELT',TRIM(reprcptr),tab) Bmelt(:,:)=tab(:,:) if (icpt.ne.3) then ! pas de lecture Hwater si icompteur = 3 call Read_ncdf_var('HWATER',TRIM(reprcptr),tab) Hwater(:,:)=tab(:,:) else Hwater(:,:)=0. end if ! lecture vitesses dans tous les cas. call Read_ncdf_var('UXBAR',TRIM(reprcptr),tab) UXBAR(:,:)=tab(:,:) call Read_ncdf_var('UYBAR',TRIM(reprcptr),tab) UYBAR(:,:)=tab(:,:) do i=1,nx-1 do j=1,ny-1 oldu(i,j) = sqrt(((uxbar(i,j)+uxbar(i+1,j))/2)**2 & + ((uybar(i,j)+uybar(i,j+1))/2.)**2) end do end do ! tableaux 3D ! ----------------------------- call Read_ncdf_var('T',TRIM(reprcptr),tab1T) T(:,:,:)=tab1T(:,:,:) call Read_ncdf_var('UX',TRIM(reprcptr),tab1) Ux(:,:,:)=tab1(:,:,:) call Read_ncdf_var('UY',TRIM(reprcptr),tab1) Uy(:,:,:)=tab1(:,:,:) if (itracebug.eq.1) call tracebug('avant lectures variables traceur') call Read_ncdf_var('XDEP',TRIM(reprcptr),tab1) !aurel neem Xdep_out(:,:,:)=tab1(:,:,:) call Read_ncdf_var('YDEP',TRIM(reprcptr),tab1) Ydep_out(:,:,:)=tab1(:,:,:) call Read_ncdf_var('TDEP',TRIM(reprcptr),tab1) Tdep_out(:,:,:)=tab1(:,:,:)+time end if ascii_nc do i=1,nx-1 do j=1,ny-1 oldu(i,j) = sqrt(((Uxbar(i,j)+Uxbar(i+1,j))/2)**2 & + ((Uybar(i,j)+Uybar(i,j+1))/2.)**2) end do end do ! moyenne de l'epaisseur et pentes do j=2,ny do i=2,nx Hmy(i,j)=(H(i,j)+H(i,j-1))/2. Hmx(i,j)=(H(i,j)+H(i-1,j))/2. Sdx(i,j)=(S(i,j)-S(i-1,j))/dx Sdy(i,j)=(S(i,j)-S(i,j-1))/dy end do end do do i=2,nx Hmx(i,1)=(H(i,1)+H(i-1,1))/2. end do do j=2,ny Hmy(1,j)=(H(1,j)+H(1,j-1))/2. end do ! vitesse a la base Uxflgz(:,:) = Ux(:,:,nz) Uyflgz(:,:) = Uy(:,:,nz) debug_3D(:,:,120) = T(:,:,nz) + 273.15 end subroutine read_recovery !> SUBROUTINE: read_no_recovery !!Pas de reprise d'un fichier compteur !> subroutine read_no_recovery use module3D_phy implicit none Real,Parameter :: Dzm=600 !< Grid Step In Mantle !Prop Thermique Real,Parameter :: Cm=1.04e8 !!!pas reprise : il faut initier les temperatures ds le socle ! temperature dans le socle lineaire avec le gradient geothermique do k=nz+1,nz+nzm do j=1,ny do i=1,nx T(i,j,k)=T(i,j,nz)-dzm*(k-nz)*ghf(i,j)/cm end do end do end do end subroutine read_no_recovery !---------------------------------------------------------- !---------------------------------------------------------- !> SUBROUTINE: out_recovery !! Sortie d'un fichier de reprise !! Attention : fichier en sortie .cptr ou fichier .nc !! used modules: - use module3D_phy !! - use netcdf !! - use io_netcdf_grisli !! - use util_recovery !! - use tracer_vars !> subroutine out_recovery(isortie) ! --------------------------------------------------------------- ! SORTIE COMPTEUR ! ! ! sortie qui stocke dans le fichier runname.cptr les variables ! S, H, B, HDOT, BDOT, et T plusieurs fois dans le programme ! (chaque nouvelle sortie ecrase la precedente) => ! on peut faire repartir le programme a partir de ce fichier, ! et donc du dernier pas de temps stocke ! ! Maintenant, il y a deux formats de fichiers possibles : ! ascii (isortie = 1) et netcdf (isortie = 2). ! --------------------------------------------------------------- use module3D_phy use netcdf use io_netcdf_grisli use util_recovery use tracer_vars ! aurel neem implicit none character(len=80) :: filin integer ncid,status real*8, dimension(:,:), pointer :: tab => null() !< tableau 2d real ecrit dans le fichier real*8, dimension(:,:,:), pointer :: tab1 => null() !< tableau 3d real real*8, dimension(:,:,:), pointer :: tab1T => null() !< tableau 3d real pour la temperature character(len=20),dimension(3) :: dimnames2d !< dimensions pour netcdf pour tableau 2d pour les noeud majeur character(len=20),dimension(4) :: dimnames3d !< pour 3d troisieme dim est nz character(len=20),dimension(4) :: dimnames3dT !< pour 3d troisieme dim est nz+nzm ! integer :: nrecs=1 !< compteur pour les enregistrements temps des variables ! integer :: idef=0 !< pour savoir si la variable a ete definie ou non integer :: isortie !< pour le choix de type de fichier en sortie if (itracebug.eq.1) call tracebug(' Entree dans routine out_recovery') call testout_recovery(filin) if (itracebug.eq.1) call tracebug(' Entree dans routine out_recovery 2') if (logic_out) then Select Case (Isortie) Case (1) ! sortie de type ascii filin=trim(filin)//'.cptr' open(num_file1,file=trim(filin)) write(num_file1,*) time_out(1), ' TIME ' write(num_file1,*) geoplace write(num_file1,*) nx*ny,nz,nzm,nx,ny write(num_file1,*) ! ecriture des tableaux 2D write(num_file1,*)' Tableaux 2D' write(num_file1,*)' Surface ' write(num_file1,*) S write(num_file1,*)' Epaisseur ' write(num_file1,*) H write(num_file1,*)' Base glace ' write(num_file1,*) B write(num_file1,*)' Socle ' write(num_file1,*) Bsoc write(num_file1,*)' Ibase ' write(num_file1,*) Ibase write(num_file1,*)' Bmelt ' write(num_file1,*) Bmelt write(num_file1,*)' Hwater ' write(num_file1,*) Hwater write(num_file1,*)' Uxbar ' write(num_file1,*) Uxbar write(num_file1,*)' Uybar ' write(num_file1,*) Uybar ! tableau 3D write(num_file1,*) write(num_file1,*) 'Temperature (3D y compris le socle)' write(num_file1,*) T(:,:,:) write(num_file1,*) write(num_file1,*) 'vitesse selon x' write(num_file1,*) Ux(:,:,:) write(num_file1,*) write(num_file1,*) 'vitesse selon y' write(num_file1,*) Uy(:,:,:) !aurel neem : traceurs dans le cptr write(num_file1,*) write(num_file1,*) 'origine de la glace en x' write(num_file1,*) xdep(:,:,:) write(num_file1,*) write(num_file1,*) 'origine de la glace en y' write(num_file1,*) ydep(:,:,:) write(num_file1,*) write(num_file1,*) 'age de la glace' tab1(:,:,:) = tdep(:,:,:) - time write(num_file1,*) tab1(:,:,:) ! note: on decide d'ecrire le temps 'relatif' close(num_file1) Case (2) ! sortie de type netcdf if (.not.associated(tab)) then allocate(tab(nx,ny),tab1(nx,ny,nz),tab1T(nx,ny,nz+nzm)) end if filin=trim(filin)//'.nc' ! creation du fichier !------------------------ if (ncdf_type.eq.32) then status = nf90_create(TRIM(filin),NF90_WRITE,ncid) ! en 32 bits else if (ncdf_type.eq.64) then status = nf90_create(trim(filin),and(nf90_write,nf90_64bit_offset),ncid) ! sur r2d2 else write(6,*)'pb de lecture de netcdf_type dans out_cptr:', ncdf_type endif status = nf90_close(ncid) ! fermeture call initfile(nx,ny,nz,nzm,dimnames2d,dimnames3d,dimnames3dT,filin) ! initialise le fichier nc call write_ncdf_var('temps','time',trim(filin),time_out,'double') ! ecrit le temps du snapshot ! variables 2D tab(:,:)= S(:,:) call write_ncdf_var('S',dimnames2d,trim(filin),tab,'double') ! ecrit S tab(:,:)= H(:,:) call write_ncdf_var('H',dimnames2d,trim(filin),tab,'double') tab(:,:)= B(:,:) call write_ncdf_var('B',dimnames2d,trim(filin),tab,'double') tab(:,:)= Bsoc(:,:) call write_ncdf_var('BSOC',dimnames2d,trim(filin),tab,'double') tab(:,:)= IBASE(:,:) call write_ncdf_var('IBASE',dimnames2d,trim(filin),tab,'double') tab(:,:)= Bmelt(:,:) call write_ncdf_var('BMELT',dimnames2d,trim(filin),tab,'double') tab(:,:)= Hwater(:,:) call write_ncdf_var('HWATER',dimnames2d,trim(filin),tab,'double') tab(:,:)= UXBAR(:,:) call write_ncdf_var('UXBAR',dimnames2d,trim(filin),tab,'double') tab(:,:)= UYBAR(:,:) call write_ncdf_var('UYBAR',dimnames2d,trim(filin),tab,'double') ! variables 3D tab1T(:,:,:)= T(:,:,:) call write_ncdf_var('T',dimnames3dT,trim(filin),tab1T,'double') ! ecrit T tab1(:,:,:)= UX(:,:,:) call write_ncdf_var('UX',dimnames3d,trim(filin),tab1,'double') tab1(:,:,:)= UY(:,:,:) call write_ncdf_var('UY',dimnames3d,trim(filin),tab1,'double') ! variables 3D pour les traceurs (Aurel NEEM) ! attention ce sont les tabelaux _out qui ont les bonnes dimensions. tab1(:,:,:)=XDEP_out(:,:,:) call write_ncdf_var('XDEP',dimnames3d,trim(filin),tab1,'double') !aurel neem tab1(:,:,:)=YDEP_out(:,:,:) call write_ncdf_var('YDEP',dimnames3d,trim(filin),tab1,'double') tab1(:,:,:)=TDEP_out(:,:,:) - time call write_ncdf_var('TDEP',dimnames3d,trim(filin),tab1,'double') End Select End if end subroutine out_recovery !> subroutine initfile !! Initialise the netcdf file !! Used module: - use netcdf !! - use io_netcdf_grisli !! @param nxx dimension along x !! @param nyy dimension along y !! @param nzz dimension along z !! @param nzmm dimension along z for T !! @param fil_sortie name of the file to initialise !! @param dimnames2d dimensions for netcdf !> subroutine initfile(nxx,nyy,nzz,nzmm,& dimnames2d,dimnames3d,dimnames3dT,file) use netcdf use io_netcdf_grisli implicit none integer,intent(in) :: nxx !< dimension along x integer,intent(in) :: nyy !< dimension along y integer,intent(in) :: nzz !< dimension along z integer,intent(in) :: nzmm !< dimension along z in socle character(len=20),dimension(3) :: dimnames2d !< dimensions pour netcdf pour tableau 2d pour les noeud majeur character(len=20),dimension(4) :: dimnames3d !< pour 3d troisieme dim est nz character(len=20),dimension(4) :: dimnames3dT !< pour 3d troisieme dim est nz+nzm character(len=*) ::file !< name of the file to init ! initialisation call write_ncdf_dim('x',trim(file),nxx) !< dimensions des tableaux call write_ncdf_dim('y',trim(file),nyy) call write_ncdf_dim('z',trim(file),nzz) call write_ncdf_dim('zm',trim(file),nzz+nzmm) call write_ncdf_dim('time',trim(file),0) dimnames2d(1)='x' dimnames2d(2)='y' dimnames2d(3)='time' dimnames3d(1)='x' dimnames3d(2)='y' dimnames3d(3)='z' dimnames3d(4)='time' dimnames3dT(1)='x' dimnames3dT(2)='y' dimnames3dT(3)='zm' dimnames3dT(4)='time' end subroutine initfile !> SUBROUTINE: symetry_cptr !! Subroutine utilisee dans les experiences axysymetriques pour repartir d'un cptr dans lequel !! tous les champs sont axysmetriques ou symetriques par rapport a un axe !! Used module: - use module3D_phy !! @param iaxe symetrie par rapport a iaxe !! @param jaxe symetrie par rapport a jaxe !> subroutine symetry_cptr(jaxe) use module3D_phy use tracer_vars implicit none integer :: jaxe integer :: jsym !symetrique par rapport à jaxe la référence est le bas. Pour l'instant seulement lui jsym=0 do j=jaxe+1,ny jsym=jsym+1 do i=1,nx S(i,j)=S(i,jaxe-jsym) H(i,j)=H(i,jaxe-jsym) B(i,j)=B(i,jaxe-jsym) Bsoc(i,j)=Bsoc(i,jaxe-jsym) Ibase(i,j)=Ibase(i,jaxe-jsym) Bmelt(i,j)=bmelt(i,jaxe-jsym) hwater(i,j)=hwater(i,jaxe-jsym) uxbar(i,j)=uxbar(i,jaxe-jsym) uybar(i,j)=-uybar(i,jaxe-jsym+1) ! attention different des autres ! a cause des staggered grids T(i,j,:)=T(i,jaxe-jsym,:) ux(i,j,:)=ux(i,jaxe-jsym,:) uy(i,j,:)=-uy(i,jaxe-jsym+1,:) ! different des autres XDEP(i,j,:)=XDEP(i,jaxe-jsym,:) YDEP(i,j,:)=YDEP(i,jaxe-jsym,:) TDEP(i,j,:)=TDEP(i,jaxe-jsym,:) end do end do return end subroutine symetry_cptr