source: trunk/SOURCES/Hudson_files/lect-hudson_mod.f90 @ 334

module lect_topo_hudson_shelf

  use module3D_phy
  use sedim_declar


  character(len=80) :: filin
  character(len=80) :: filin2

contains
 
  subroutine input_topo


    write(6,*)' Heino: entree dans routine input_topo'
!====================================== La reponse est 42 ===========
    write(42,*)
    write(42,*)' Fichiers en entree'
    write(42,*)'----------------------'
!==================================================================== 

    nxx=nx
    nyy=ny
    print*,'nx,ny',nx,ny

!jalv: nouvelle facon de determiner la topo: lecture du mask sediments
!      et mask topographie independement.


    filin='mask-sedim-hudson.dat'

    filin2='topo-init-hudson.dat'

    filin = TRIM(DIRNAMEINP)//TRIM(filin)

    filin2 = TRIM(DIRNAMEINP)//TRIM(filin2)



!............TOPOGRAPHIE....HUDSON...................................

! jalv: ouverture du  masque topo :
    open(num_coupe,file=filin2)
    do i=1,6
       read(num_coupe,*) ! 6 lignes de commentaires
    end do
! lecture proprement dite
    do j=ny, 1, -1
       read(num_coupe,'(151i1)') (mk(i,j), i=1,nx)
    end do
    close(num_coupe)


    do j=1,NY
       do i=1,NX
         IF (MK(i,j).EQ.9) THEN ! Surface a 400 metres


             S(i,j)=400.
             S0(i,j)=400.
             BSOC(i,j)=400.
             B(i,j)=400.
         ENDIF


         IF (MK(i,j).EQ.8) THEN! Surface a 100 mettres 
             S(i,j)=100.
             S0(i,j)=100.
             BSOC(i,j)=100.
             B(i,j)=100.
         ENDIF

         IF (MK(i,j).EQ.7) THEN! ocean  peu profond (200 m)
             S(i,j)=0.
             S0(i,j)=0.
             BSOC(i,j)=-200.
             B(i,j)=0.
         ENDIF

         IF (MK(i,j).EQ.6) THEN! ocean interm (600 m)
             S(i,j)=0.
             S0(i,j)=0.
             BSOC(i,j)=-600.
             B(i,j)=0.
         ENDIF
         IF (MK(i,j).EQ.5) THEN! ocean (1200 m)
             S(i,j)=0.
             S0(i,j)=0.
             BSOC(i,j)=-1200.
             B(i,j)=0.
         ENDIF

         IF (MK(i,j).EQ.4) THEN! ocean  (1700 m)
             S(i,j)=0.
             S0(i,j)=0.
             BSOC(i,j)=-1800.
             B(i,j)=0.
         ENDIF

         IF (MK(i,j).EQ.3) THEN! ocean  (2200 m)
             S(i,j)=0.
             S0(i,j)=0.
             BSOC(i,j)=-2400.
             B(i,j)=0.
         ENDIF

         IF (MK(i,j).EQ.1) THEN! ocean tres profond (2800 m)
             S(i,j)=0.
             S0(i,j)=0.
             BSOC(i,j)=-2800.
             B(i,j)=0.
         ENDIF

       enddo
    enddo

!......FINI.....TOPO...............................................



!.............MASK...SEDIMENTS....................................


!jalv:  ouverture du  masque sediments :
    open(num_coupe,file=filin)
    do i=1,6
       read(num_coupe,*) ! 6 lignes de commentaires
    end do
! lecture proprement dite
    do j=ny, 1, -1
       read(num_coupe,'(151i1)') (mk(i,j), i=1,nx)
    end do
    close(num_coupe)

! initialisation du mask de sediment mksedim :
    mksedim(:,:)=mk(:,:)
! a la frontiere le demi noeud est sediment.
! calcul du masque sur les demi mailles : frontiere -> sedim
    mkxsedim(:,:)=max(mksedim(:,:),eoshift(mksedim(:,:),shift=-1,boundary=0,dim=1))
    mkysedim(:,:)=max(mksedim(:,:),eoshift(mksedim(:,:),shift=-1,boundary=0,dim=2))
!    write(42,*) 'sur les bords zone sediment, masque exterieur '


    where (mk(:,:)>1) mk(:,:)=0         ! le masque mk est terre (0) ou ocean (1)


    where (BSOC(:,:).le.0) mk(:,:)=1       ! jalv: pour qu'il voit l'ocean lors de l'initialisation 



!..................................................................
!.........FINI.....MASQUE....SEDIMENTS.............................



! masque mko vrai partout pour autoriser la presence de glace dans icethick
! Attention MK0 defini les points de grille ou on autrorise la presence de glace
    MK0(:,:)=0      
    do I=2,NX-1
       do J=2,NY-1
          mk0(I,J)=1
       end do
    end do
!call init_sliding

    xlong(:,:)=0.
    ylat(:,:)=0.
    
    xmin=0.
    ymin=0.
    xmax=(nx-1)*dx/1000.
    ymax=(ny-1)*dx/1000.

!write(num_rep_42,*) 'domaine geographique en km, 0 en bas gauche',int(xmax),'x',int(ymax)
    write(42,*) 'domaine geographique en km : xmin : ',int(xmin),' ymin : ',int(ymin), &
                                                 ' xmax : ',int(xmax),' ymax : ',int(ymax)

! flux geothermique 42 mW/m2
    ghf(:,:)=42.

    write(42,*) 'flux geothermique uniforme (mW/m2) :',ghf(1,1)
    ghf(:,:)=-SECYEAR/1000.*ghf(:,:)



  end subroutine input_topo


end module lect_topo_hudson_shelf