module lect_topo_snowball use module3D_phy character(len=100) :: topo_ref character(len=100) :: topo_dep character(len=100) :: grid_topo real :: toto real :: correc_bathy contains subroutine input_topo namelist/topo_file/topo_ref,topo_dep,grid_topo 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 ! nom du bloc ' write(num_rep_42,*) write(num_rep_42,*) 'topo_ref = ', topo_ref write(num_rep_42,*) 'topo_dep = ', topo_dep write(num_rep_42,*) 'grid_topo = ', grid_topo 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,*) ! file1=TRIM(DIRNAMEINP)//'topo-toanor.g40' ! topo-maasnor.g40' ! topo de depart ! file2=TRIM(DIRNAMEINP)//'topo-toanor.g40' !'topo-maasnor.g40' ! topo de reference (isostasie) ! write(42,*) 'topo de depart', file1 ! write(42,*) 'topo reference', file2 ! lecture de la topo de reference ! --------------------------- open (20,file=TRIM(DIRNAMEINP)//trim(topo_ref)) read(20,'(A80)') TITRE read(20,*) NI,NJ,NXX,NYY,STEP read(20,*) do J=1,ny do I=1,nx read (20,*) toto,toto,BSOC0(I,J) H0(i,j)=0. S0(i,j)=max(BSOC0(i,j),0.) end do end do close(20) ! pour ne pas avoir de bathymetrie < -5000m ! where (BSOC0(:,:).LT.-5000) BSOC0(:,:)=-5000. ! lecture de la topo de depart ! --------------------------- open (20,file=TRIM(DIRNAMEINP)//trim(topo_dep)) ! 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,*) toto,toto,BSOC(I,J) H(i,j)=0. S(i,j)=max(BSOC(i,j),0.) B(i,j)=S(i,j)-H(i,j) end do end do close(20) ! correction de la bathymetrie : where (BSOC(:,:).LT.0.) BSOC(:,:)=-5000. correc_bathy=0. do j=2,ny-1 do i=2,nx-1 if (BSOC(i,j).LT.0.) THEN if (BSOC(i+1,j).gt.correc_bathy.or.BSOC(i,j+1).gt.correc_bathy.or.BSOC(i-1,j).gt.correc_bathy.or.bsoc(i,j-1).gt.correc_bathy.or.bsoc(i+1,j+1).gt.correc_bathy.or.bsoc(i-1,j-1).gt.correc_bathy) then BSOC(i,j)=-50. endif endif enddo enddo correc_bathy=-50. do j=2,ny-1 do i=2,nx-1 if (BSOC(i,j).LT.correc_bathy) THEN if (BSOC(i+1,j).eq.correc_bathy.or.BSOC(i,j+1).eq.correc_bathy.or.BSOC(i-1,j).eq.correc_bathy.or.bsoc(i,j-1).eq.correc_bathy.or.bsoc(i+1,j+1).eq.correc_bathy.or.bsoc(i-1,j-1).eq.correc_bathy) then BSOC(i,j)=-200. endif endif enddo enddo correc_bathy=-200. do j=2,ny-1 do i=2,nx-1 if (BSOC(i,j).LT.correc_bathy) THEN if (BSOC(i+1,j).eq.correc_bathy.or.BSOC(i,j+1).eq.correc_bathy.or.BSOC(i-1,j).eq.correc_bathy.or.bsoc(i,j-1).eq.correc_bathy.or.bsoc(i+1,j+1).eq.correc_bathy.or.bsoc(i-1,j-1).eq.correc_bathy) then BSOC(i,j)=-500. endif endif enddo enddo correc_bathy=-500. do j=2,ny-1 do i=2,nx-1 if (BSOC(i,j).LT.correc_bathy) THEN if (BSOC(i+1,j).eq.correc_bathy.or.BSOC(i,j+1).eq.correc_bathy.or.BSOC(i-1,j).eq.correc_bathy.or.bsoc(i,j-1).eq.correc_bathy.or.bsoc(i+1,j+1).eq.correc_bathy.or.bsoc(i-1,j-1).eq.correc_bathy) then BSOC(i,j)=-750. endif endif enddo enddo correc_bathy=-750. do j=2,ny-1 do i=2,nx-1 if (BSOC(i,j).LT.correc_bathy) THEN if (BSOC(i+1,j).eq.correc_bathy.or.BSOC(i,j+1).eq.correc_bathy.or.BSOC(i-1,j).eq.correc_bathy.or.bsoc(i,j-1).eq.correc_bathy.or.bsoc(i+1,j+1).eq.correc_bathy.or.bsoc(i-1,j-1).eq.correc_bathy) then BSOC(i,j)=-1000. endif endif enddo enddo correc_bathy=-1000. do j=2,ny-1 do i=2,nx-1 if (BSOC(i,j).LT.correc_bathy) THEN if (BSOC(i+1,j).eq.correc_bathy.or.BSOC(i,j+1).eq.correc_bathy.or.BSOC(i-1,j).eq.correc_bathy.or.bsoc(i,j-1).eq.correc_bathy.or.bsoc(i+1,j+1).eq.correc_bathy.or.bsoc(i-1,j-1).eq.correc_bathy) then BSOC(i,j)=-1500. endif endif enddo enddo correc_bathy=-1500. do j=2,ny-1 do i=2,nx-1 if (BSOC(i,j).LT.correc_bathy) THEN if (BSOC(i+1,j).eq.correc_bathy.or.BSOC(i,j+1).eq.correc_bathy.or.BSOC(i-1,j).eq.correc_bathy.or.bsoc(i,j-1).eq.correc_bathy.or.bsoc(i+1,j+1).eq.correc_bathy.or.bsoc(i-1,j-1).eq.correc_bathy) then BSOC(i,j)=-2000. endif endif enddo enddo correc_bathy=-2000. do j=2,ny-1 do i=2,nx-1 if (BSOC(i,j).LT.correc_bathy) THEN if (BSOC(i+1,j).eq.correc_bathy.or.BSOC(i,j+1).eq.correc_bathy.or.BSOC(i-1,j).eq.correc_bathy.or.bsoc(i,j-1).eq.correc_bathy.or.bsoc(i+1,j+1).eq.correc_bathy.or.bsoc(i-1,j-1).eq.correc_bathy) then BSOC(i,j)=-3000. endif endif enddo enddo correc_bathy=-3000. do j=2,ny-1 do i=2,nx-1 if (BSOC(i,j).LT.correc_bathy) THEN if (BSOC(i+1,j).eq.correc_bathy.or.BSOC(i,j+1).eq.correc_bathy.or.BSOC(i-1,j).eq.correc_bathy.or.bsoc(i,j-1).eq.correc_bathy.or.bsoc(i+1,j+1).eq.correc_bathy.or.bsoc(i-1,j-1).eq.correc_bathy) then BSOC(i,j)=-4000. endif endif enddo enddo ! pour ne pas avoir de bathymetrie aberante : where (BSOC(:,:).LT.-5000) BSOC(:,:)=-5000. ! bathymetrie de reference : BSOC0(:,:)=BSOC(:,:) ! lecture des coordonnées geographiques ! les coordonnees sont calculees en °dec avec GMT, ! les longitudes sont comprises entre 0 et 360 open(unit=2004,file=trim(dirnameinp)//trim(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) write(42,*) 'fichier grille: ', filin ! passage en longitude -180 +180 : where (XLONG(:,:).GT.180) XLONG(:,:)=XLONG(:,:)-360. xmin=xcc(1,1) ymin=ycc(1,1) xmax=xcc(nx,ny) ymax=ycc(nx,ny) ! flux geothermique fixe : ghf(:,:)=-SECYEAR/1000.*42. ! pour passer les flux des mW/m2 au J/m2/an ! ghf(:,:)=-SECYEAR/1000.*ghf(:,:) ! write(42,*) 'flux geothermique fixe : 55 mW/m2' ! ghf(:,:)=-SECYEAR/1000.*55. !B6norcg2 ! print*,'lect topo' ! print*,'shb0',S0(300,150),H0(300,150),BSOC0(300,150) ! Initialisation du Masque !------------------------------------------------ ! pour l'Hemisphere Nord mko vrai partout (version 2006) MK0(:,:)=1 !------------------------------------------------ end subroutine input_topo end module lect_topo_snowball