!> \file lect-anteis_mod.f90 !!Module pour la lecture de la topographie !< !> \namespace lect_topo_anteis !! Module pour la lecture de la topographie !! \author ... !! \date ... !! @note Used module !! @note - use module3D_phy !! @note - use param_phy_mod !< module lect_topo_anteis use module3D_phy use interface_input use io_netcdf_grisli implicit none character(len=100) :: topo_dep ! Topo de départ character(len=100) :: topo_ref ! Topo de référence character(len=100) :: grid_topo ! fichier grille character(len=100) :: ghf_fich ! fichier grille character(len=100) :: filin character(len=100) :: file_ncdf !< fichier netcdf issue des fichiers .dat ! character(len=100) :: file1 ! character(len=100) :: file2 real, dimension(nx,ny,5) :: bidon ! pour l'appel a courbure real :: sealevel0 contains subroutine input_topo integer :: ios ! pour les lectures ncdf real*8, dimension(:,:), pointer :: tab !< tableau 2d real pointer namelist/topo_file/topo_ref,topo_dep,grid_topo,ghf_fich rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,topo_file) ! formats pour les ecritures dans 42 428 format(A) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '&topo_file ! input_topo ' write(num_rep_42,'(A,A)') 'topo_ref = ', topo_ref write(num_rep_42,'(A,A)') 'topo_dep = ', topo_dep write(num_rep_42,'(A,A)') 'grid_topo =', grid_topo write(num_rep_42,'(A,A)') 'ghf_fich = ', ghf_fich write(num_rep_42,*)'/' write(num_rep_42,428) '! topo_ref= topo ref isostasie' write(num_rep_42,428) '! topo_dep= topo de depart' write(num_rep_42,428) '! grid_topo : fichier i,j,x,y,lon,lat' write(num_rep_42,428) '! ghf_fich : fichier flux geothermique' write(num_rep_42,*) topo_ref=trim(dirnameinp)//trim(topo_ref) topo_dep=trim(dirnameinp)//trim(topo_dep) grid_topo=trim(dirnameinp)//trim(grid_topo) ghf_fich=trim(dirnameinp)//trim(ghf_fich) ! lecture adaptee aux fichiers intercomparaison EISMINT nxx=nx nyy=ny !!$! lecture de la topo actuelle !!$! --------------------------- !!$ open (20,file=TRIM(DIRNAMEINP)//file2,status='old') !!$ !!$ read(20,'(A80)') TITRE !!$ read(20,*) NI,NJ,NXX,NYY,STEP !!$ read(20,*) !!$ do J=1,ny !!$ do I=1,nx !!$ read (20,*) S0(I,J),H0(I,J),BSOC0(I,J) !!$ S0(i,j)=max(S0(i,j),0.) !!$ end do !!$ end do !!$ close(20) !!$ !!$ !!$! lecture de la topo de depart !!$! --------------------------- !!$ open (20,file=TRIM(DIRNAMEINP)//file1,status='old') !!$! open (20,file='../INPUT-DATA/hemin.g50') !!$ read(20,'(A80)') TITRE !!$ read(20,*) NI,NJ,NXX,NYY,STEP !!$ read(20,*) !!$ do J=1,ny !!$ do I=1,nx !!$ read (20,*) S(I,J),H(I,J),BSOC(I,J) !!$ end do !!$ end do !!$ close(20) ! lecture de la topo de référence ! ------------------------------- ! Cette topo sert a calculer le socle de reference pour l'isostasie ! voir init_iso et a avoir une surface de reference pour les temperatures ! lecture adaptee aux fichiers ZBL.dat ou netcdf ou grd ! call lect_input(1,'Bsoc',1,Bsoc0,topo_ref,file_ncdf) ! socle ! call lect_input(1,'S',1,S0,topo_ref,file_ncdf) ! surface ! call lect_input(1,'H',1,H0,topo_ref,file_ncdf) ! epaisseur ! lecture pour eviter plantage avec compile -O0 call Read_Ncdf_var('Bsoc',topo_ref,tab) Bsoc0(:,:) = tab(:,:) call Read_Ncdf_var('S',topo_ref,tab) S0(:,:) = tab(:,:) call Read_Ncdf_var('H',topo_ref,tab) H0(:,:) = tab(:,:) !~ !cdc correction de 3 pts qui posent probleme (pente tres forte): !~ S0(38,53)=1400. !~ Bsoc0(38,53)=-1000. !~ H0(38,53)=S0(38,53)-Bsoc0(38,53) !~ !~ S0(35,53)=1300. !~ Bsoc0(35,53)=-1000. !~ H0(35,53)=S0(35,53)-Bsoc0(35,53) !~ !~ S0(35,56)=1200. !~ Bsoc0(35,56)=-1000. !~ H0(35,56)=S0(35,56)-Bsoc0(35,56) !cdc correction point pole sud : ! S0(71,71)=(S0(71,70)+S0(71,72)+S0(70,71)+S0(72,71))/4. ! Bsoc0(71,71)=(Bsoc0(71,70)+Bsoc0(71,72)+Bsoc0(70,71)+Bsoc0(72,71))/4. ! where (S0(:,:).GT.0) ! H0(:,:)=S0(:,:)-BSOC0(:,:) ! elsewhere ! H0(:,:)=1. ! endwhere where (BSOC(:,:).LT.-9999.) BSOC(:,:)=-9999. endwhere sealevel0=0. ! voir a passer dans le fichier parametre S0(:,:)=max(S0(:,:),sealevel0) ! pour etre au niveau des mers : ATTENTION si SEALEV <0 mk0(:,:) = 1 ! mk0=0 pour les zones interdites ! lecture de la topo de depart ! --------------------------- ! lecture adaptee aux fichiers ZBL.dat ou netcdf ou grd ! call lect_input(1,'Bsoc',1,Bsoc,topo_dep,file_ncdf) ! socle ! call lect_input(1,'S',1,S,topo_dep,file_ncdf) ! surface ! call lect_input(1,'H',1,H,topo_dep,file_ncdf) ! epaisseur ! lecture pour eviter plantage avec compile -O0 call Read_Ncdf_var('Bsoc',topo_dep,tab) Bsoc(:,:) = tab(:,:) call Read_Ncdf_var('S',topo_dep,tab) S(:,:) = tab(:,:) call Read_Ncdf_var('H',topo_dep,tab) H(:,:) = tab(:,:) !~ !cdc correction de 3 pts qui posent probleme (pente tres forte): !~ S(39,54)=1400. !~ Bsoc(39,54)=-1000. !~ H(39,54)=S(39,54)-Bsoc(39,54) !~ !~ S(36,54)=1300. !~ Bsoc(36,54)=-1000. !~ H(36,54)=S(36,54)-Bsoc(36,54) !~ !~ S(36,57)=1200. !~ Bsoc(36,57)=-1000. !~ H(36,57)=S(36,57)-Bsoc(36,57) ! S(71,71)=(S(71,70)+S(71,72)+S(70,71)+S(72,71))/4. ! Bsoc(71,71)=(Bsoc(71,70)+Bsoc(71,72)+Bsoc(70,71)+Bsoc(72,71))/4. ! where (S(:,:).GT.0) ! H(:,:)=S(:,:)-BSOC(:,:) ! elsewhere ! H(:,:)=1. ! endwhere S(:,:)=max(S(:,:),0.) ! pour etre au niveau des mers : ATTENTION si SEALEV <0 H(:,:)=max(H(:,:),1.) ! pour avoir au moins 1 m !!$ !!$ ! socle !!$ filin='bedelev-2000-40km.dat' !!$ call lect_eis(nx,ny,BSOC,filin,DIRNAMEINP) !!$ write(num_rep_42,*) 'fichier BSOC (socle) : ', filin !!$ Bsoc0(:,:)=Bsoc(:,:) !!$ ! surface !!$ filin='surface-2000-40km.dat' !!$ call lect_eis(nx,ny,S,filin,DIRNAMEINP) !!$ write(num_rep_42,*) 'fichier surface : ', filin !!$ !!$ ! epaisseur !!$ filin='icethic-2000-40km.dat' !!$ call lect_eis(nx,ny,H,filin,DIRNAMEINP) !!$ write(num_rep_42,*) 'fichier epaisseur : ', filin !!$ H0(:,:)=H(:,:) !!$ !!$ !!$ ! masque !!$ filin='maskHUY40km.dat' !!$ call lect_ieis(nx,ny,MK0,filin,DIRNAMEINP) !!$ write(num_rep_42,*) 'fichier masque : ', filin !!$ !!$ ! enlever les epaisseurs fictives de Philippe !!$ do I=1,NX !!$ do J=1,NY !!$ ! attention les valeurs de masque de Philippe sont differentes !!$ if ((MK0(I,J).gt.1).and.(H(i,j).le.218.00)) then ! mer libre !!$ H(i,j)=1. !!$ S(I,J)=H(i,j)-(RO/ROW)*H(i,j) !!$ endif !!$ end do !!$ end do ! S0(:,:)=S(:,:) ! pour l'Antarctique masque mko vrai partout (version 2006) MK0(:,:)=1 ! determination des flot do J=1,NY do I=1,NX if (((BSOC(I,J)+H(I,J)*RO/ROW -SEALEVEL).LT.0.).and.(H(I,J).gt.1.E-3)) then FLOT(I,J)=.TRUE. else FLOT(I,J)=.FALSE. endif enddo enddo ! calcul des courbures du socle ! call courbure(nx,ny,dx,Bsoc,bidon(:,:,1),bidon(:,:,2),bidon(:,:,3), & ! bidon(:,:,4),socle_cry,bidon(:,:,5)) ! socle_cry(:,:)=socle_cry(:,:)*dx*dx ! lecture des coordonnées geographiques ! filin=TRIM(DIRNAMEINP)//'coord-Ant-40km.dat' ! les coordonnees sont calculees en °dec avec GMT, ! les longitudes sont comprises entre -180 et +180 (negative a l'Ouest de ! Greenwich et positive a l'Est) open(unit=2004,file=grid_topo,iostat=ios) do k=1,nx*ny read(2004,*) i,j,xcc(i,j),ycc(i,j),Xlong(i,j),Ylat(i,j) enddo close(2004) xmin=xcc(1,1)/1000. ymin=ycc(1,1)/1000. xmax=xcc(nx,ny)/1000. ymax=ycc(nx,ny)/1000. !!$ ! lecture du fichier de reference pour le calcul du niveau des mers : etat actuel !!$ open(88,file=TRIM(DIRNAMEINP)//'grzone56-k000.pl') !!$ !!$ read(88,*) !!$ read(88,*) !!$ read(88,*) !!$ do j=1,ny !!$ do i=1,nx !!$ read(88,*) S_sealev(i,j),H_sealev(i,j),B_sealev(i,j),M_sealev(i,j) !!$ enddo !!$ enddo !!$ close(88) !!$ !!$ write(num_rep_42,*) 'fichier reference pour le niveau des mers : ', & !!$ TRIM(DIRNAMEINP)//'grzone56-k000.pl' ! lecture du flux geothermique de Shapiro open(88,file=ghf_fich) ! write(num_rep_42,*) 'flux geothermique Shapiro : ', TRIM(DIRNAMEINP)//'ijphi-40km-ant.dat' ! do k=1,nx*ny ! read(88,*) i,j,ghf(i,j) ! end do ! close(88) ! call lect_input(1,'ghf',1,ghf,ghf_fich,file_ncdf) ! pour eviter plantage -O0 call Read_Ncdf_var('ghf',ghf_fich,tab) ghf(:,:) = tab(:,:) ! pour passer les flux des mW/m2 au J/m2/an ghf(:,:)=-SECYEAR/1000.*ghf(:,:) !------------------------------------------------ ! mko vrai partout (version 2006) MK0(:,:)=1 !------------------------------------------------ end subroutine input_topo end module lect_topo_anteis