module dragging_hwat_sedim ! Defini les zones de stream avec : ! * un critere sur la hauteur d'eau ! * un critere sur la courbure du socle (si negatif, vallees) ! * un critere sur l'epaisseur de sediments (carte de Gabi & Laske) use module3d_phy implicit none logical,dimension(nx,ny) :: fleuvemx logical,dimension(nx,ny) :: fleuvemy logical,dimension(nx,ny) :: fleuve logical,dimension(nx,ny) :: fleuve_sedim logical,dimension(nx,ny) :: fleuve_sedimx logical,dimension(nx,ny) :: fleuve_sedimy logical,dimension(nx,ny) :: cote real :: valmax integer :: imax,jmax integer :: i_moins1,i_plus1,j_moins1,j_plus1 integer :: lmax=20 integer :: idep,jdep,iloc,jloc real :: tostick ! pour la glace posee real :: tob_ile ! pour les iles real :: cry_lim=50. ! courbure limite pour le suivi des fleuves real,dimension(nx,ny) :: h_sedim ! epaisseur de sediment real,dimension(nx,ny) :: neff ! pression effective noeuds majeurs real :: seuil_hwater ! seuil sur hwater pour avoir du glissement en zone sediment real :: seuil_sedim ! seuil sur h_sedim pour avoir du glissement en zone sediment real :: seuil_neff ! seuil sur la pression effective pour avoir glissement en zone sediment real :: coef_sedim ! coef frottement zones sediments (1) real :: coef_gz ! coef frottement zones stream std (10) real :: coef_ile ! coef frottement zones iles (0.1) real :: toto real :: betacotemax=100 ! valeur max pour betamx et betamy sur les cote real :: betacotemin=4 real :: betasedmax=400 ! 400 dans c-HadCM3-evol-basal-melt-sed : cHadCM3s et 250 dans c-HadCM3-evol-basal-melt-sed-betasedimin6 real :: betasedmin=10 ! 10 dans c-HadCM3-evol-basal-melt-sed : cHadCM3s et 6 dans c-HadCM3-evol-basal-melt-sed-betasedimin6 real :: betastreamax=500 real :: betastreamin=10 real :: betailemax=250 real :: betailemin=6 real, dimension(nx,ny) :: Vcol_x !< uniquement pour compatibilite avec spinup cat real, dimension(nx,ny) :: Vcol_y !< uniquement pour compatibilite avec spinup cat real, dimension(nx,ny) :: Vsl_x !< uniquement pour compatibilite avec spinup cat real, dimension(nx,ny) :: Vsl_y !< uniquement pour compatibilite avec spinup cat logical :: corr_def = .false. !< for deformation correction, pour compatibilite beta character(len=150) :: file_sedim ! fichier epaisseur sediment contains !------------------------------------------------------------------------------- subroutine init_dragging ! Cette routine fait l'initialisation du dragging. implicit none namelist/drag_hwat_sedim/file_sedim,hwatstream,cf,betamax,toblim,seuil_hwater,seuil_sedim,seuil_neff,coef_sedim,coef_gz,coef_ile if (itracebug.eq.1) write(num_tracebug,*)' dragging avec hwatermax' iter_beta=0 ! formats pour les ecritures dans 42 428 format(A) ! lecture des parametres du run block draghwat !-------------------------------------------------------------------- rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,drag_hwat_sedim) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '&drag_hwat_sedim ! nom du bloc hwat_sedim ' write(num_rep_42,*) write(num_rep_42,*) 'file_sedim = ', trim(file_sedim) write(num_rep_42,*) 'hwatstream = ', hwatstream write(num_rep_42,*) 'cf = ', cf write(num_rep_42,*) 'betamax = ', betamax write(num_rep_42,*) 'toblim = ', toblim write(num_rep_42,*) 'seuil_hwater = ', seuil_hwater write(num_rep_42,*) 'seuil_sedim = ', seuil_sedim write(num_rep_42,*) 'seuil_neff = ', seuil_neff write(num_rep_42,*) 'coef_sedim = ', coef_sedim write(num_rep_42,*) 'coef_gz = ', coef_gz write(num_rep_42,*) 'coef_ile = ', coef_ile write(num_rep_42,*) '/' write(num_rep_42,428) '! file_sedim : fichier carte epaisseur sediment' write(num_rep_42,428) '! hwatstream (m) : critere de passage en stream en partant de la cote si hwater > hwatstream ' write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff' write(num_rep_42,428) '! seulement pour les points cotiers' write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams ' write(num_rep_42,428) '! toblim : (Pa) pour les iles ' write(num_rep_42,428) '! seuil_hwater (m) : seuil hwater pour avoir glissement sur zone sediment' write(num_rep_42,428) '! seuil_sedim (m) : seuil epaisseur sediment pour avoir glissement sur zone sediment' write(num_rep_42,428) '! seuil_neff (Pa) seuil sur la pression effective pour avoir glissement en zone sediment' write(num_rep_42,428) '! coef_sedim : coef frottement zones stream sediments' write(num_rep_42,428) '! coef_gz : coef frottement zones stream std' write(num_rep_42,428) '! coef_ile : coef frottement zones iles' write(num_rep_42,*) tostick=1.e5 ! valeurs pour les points non flgzmx tob_ile=betamax/2. moteurmax=toblim !------------------------------------------------------------------- ! masque stream mstream_mx(:,:)=1 mstream_my(:,:)=1 ! coefficient permettant de modifier le basal drag. drag_mx(:,:)=1. drag_my(:,:)=1. where (socle_cry(:,:).lt.cry_lim) debug_3D(:,:,25)=1 elsewhere debug_3D(:,:,25)=0 endwhere ! lecture du fichier sediment : !open(20,file=TRIM(DIRNAMEINP)//trim(file_sedim)) ! 'sediment_ij_hemin40.dat') !do j=1,ny ! do i=1,nx ! read(20,*),toto,toto,h_sedim(i,j) ! enddo !enddo !h_sedim(:,:)=max(h_sedim(:,:),0.) ! pas d'epaisseurs de sediments negatives h_sedim(:,:)=10. return end subroutine init_dragging !________________________________________________________________________________ !------------------------------------------------------------------------- subroutine dragging ! defini la localisation des streams et le frottement basal ! les iles n'ont pas de condition neff mais ont la condition toblim ! (ce bloc etait avant dans flottab) do j=2,ny do i=2,nx ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim) ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim) end do end do ! le calcul des fleuves se fait sur les mailles majeures en partant de la cote ! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line ! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max ! tant que ce point est superieur a hwatermax. ! Attention : ce systeme ne permet pas d'avoir des fleuves convergents ! et les fleuves n'ont une épaisseur que de 1 noeud. fleuvemx(:,:)=.false. fleuvemy(:,:)=.false. fleuve(:,:)=.false. cote(:,:)=.false. ! detection des cotes do j=2,ny-1 do i=2,nx-1 if ((.not.flot(i,j)).and. & ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then cote(i,j)=.true. endif end do end do ! calcul de neff (pression effective sur noeuds majeurs) do j=1,ny-1 do i=1,nx-1 neff(i,j)=(neffmx(i,j)+neffmx(i,j+1)+neffmy(i,j)+neffmy(i+1,j))/4 enddo enddo ! test sur hauteur d'eau et epaisseur sediment : !where (h_sedim(:,:).GT.seuil_sedim.and.hwater(:,:).GT.seuil_hwater.and.(.not.flot(:,:))) fleuve(:,:)=.true. where (h_sedim(:,:).GT.seuil_sedim.and.neff(:,:).LE.seuil_neff.and.hwater(:,:).GT.seuil_hwater.and.(.not.flot(:,:))) fleuve(:,:)=.true. fleuve_sedim(:,:)=fleuve(:,:) ! points uniquement fleuves sediment fleuve_maj: do j=2,ny-1 ifleuve: do i=2,nx-1 cote_detect : if (cote(i,j)) then idep=i jdep=j if (socle_cry(i,j).lt.0.) then ! dans une vallee fleuve(i,j)=.true. ! write(166,*) ! write(166,*)'cote',i,j,hwater(i,j),hwatstream else cote(i,j)=.false. cycle ifleuve endif suit : do l=1,lmax ! debut de la boucle de suivi, lmax longueur maxi des fleuves i_moins1=max(idep-1,2) j_moins1=max(jdep-1,1) i_plus1=min(idep+1,nx) j_plus1=min(jdep+1,ny) !!$! recherche du max en suivant la hauteur d'eau maxi chez les voisins !!$! * en excluant les points flottants !!$! * et ceux qui sont deja tagges fleuves !version 18 juillet 2007 !!$ valmax=0. !!$ !!$ do jloc=j_moins1,j_plus1 !!$ do iloc=i_moins1,i_plus1 !!$ !!$ if ((hwater(iloc,jloc).gt.valmax) & !!$ .and.(.not.flot(iloc,jloc)) & !!$ .and.(.not.fleuve(iloc,jloc))) then !!$ imax=iloc !!$ jmax=jloc !!$ valmax=hwater(iloc,jloc) !!$ endif !!$ end do !!$ end do ! recherche du max en suivant le socle le plus profond ! * en excluant les points flottants ! * et ceux qui sont deja tagges fleuves valmax=1000. do jloc=j_moins1,j_plus1 do iloc=i_moins1,i_plus1 if ((B(iloc,jloc).lt.valmax) & .and.(.not.flot(iloc,jloc)) & .and.(.not.fleuve(iloc,jloc)).and.(socle_cry(iloc,jloc).lt.cry_lim)) then imax=iloc jmax=jloc valmax=B(iloc,jloc) endif end do end do if ((hwater(imax,jmax).gt.hwatstream).and.(socle_cry(i,j).lt.cry_lim)) then fleuve(imax,jmax)=.true. idep=imax jdep=jmax ! write(166,*) imax,jmax,hwater(imax,jmax) else ! fleuve(imax,jmax)=.false. exit suit end if end do suit end if cote_detect end do ifleuve end do fleuve_maj ! detection de problemes !write(166,*) !write(166,*) 'detection de pb dans fleuve time=',time !do j=2,ny ! do i=2,ny ! if (fleuve(i,j).and.(socle_cry(i,j).gt.cry_lim)) then ! write(166,*) 'pb fleuve',time,i,j,hwater(i,j),socle_cry(i,j) ! endif ! end do !end do do j=1,ny-1 do i=1,nx-1 if (fleuve(i,j)) then fleuvemx(i,j)=.true. fleuvemx(i+1,j)=.true. fleuvemy(i,j)=.true. fleuvemy(i,j+1)=.true. end if if (fleuve_sedim(i,j)) then fleuve_sedimx(i,j)=.true. fleuve_sedimx(i+1,j)=.true. fleuve_sedimy(i,j)=.true. fleuve_sedimy(i,j+1)=.true. endif end do end do ! pas de fleuve dans les endroits trop en pente !fleuvemx(:,:)=fleuvemx(:,:).and.(abs(rog*Hmx(:,:)*sdx(:,:)).lt.toblim) !fleuvemy(:,:)=fleuvemy(:,:).and.(abs(rog*Hmy(:,:)*sdy(:,:)).lt.toblim) !fleuvemx(:,:)=fleuvemx(:,:).and.(abs(rog*Hmx(:,:)*sdx(:,:)).lt.toblim) !fleuvemy(:,:)=fleuvemy(:,:).and.(abs(rog*Hmy(:,:)*sdy(:,:)).lt.toblim) do j=1,ny do i=1,nx if ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) gzmx(i,j)=.true. ! calcul du frottement basal (ce bloc etait avant dans neffect) if (cotemx(i,j)) then ! point cotier betamx(i,j)=cf*neffmx(i,j)*0.1 ! cat betamx(i,j)=min(betamx(i,j),tobmax) ! version Jorge ! point cotier: frott maxi=40, mini=4 ! betamx(i,j)=min(betamx(i,j),betacotemax) ! betamx(i,j)=max(betamx(i,j),betacotemin) ! version Jorge stream en zone sediment : else if (fleuve_sedimx(i,j).and.gzmx(i,j).and.(.not.cotemx(i,j))) then ! stream en zone sediment ! stream sediment: frott maxi=60, mini=6 !cdc betamx(i,j)=cf*neffmx(i,j)*3 betamx(i,j)=cf*neffmx(i,j)*coef_sedim ! betamx(i,j)=min(betamx(i,j),betasedmax) ! betamx(i,j)=max(betamx(i,j),betasedmin) else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then ! stream normal ! cat betamx(i,j)=tobmax ! cat betamx(i,j)=cf*neffmx(i,j)*20. ! attention a cette modif juillet 2007 ! cat betamx(i,j)=min(betamx(i,j),tobmax) ! cat betamx(i,j)=max(betamx(i,j),10.) ! version Jorge ! stream normal: frott maxi=80, mini=8 betamx(i,j)=cf*neffmx(i,j)*coef_gz ! betamx(i,j)=min(betamx(i,j),betastreamax) ! betamx(i,j)=max(betamx(i,j),betastreamin) ! cat if (cf*neffmx(i,j).gt.1500.) then ! cat gzmx(i,j)=.false. ! cat betamx(i,j)=tostick ! cat endif else if (ilemx(i,j)) then betamx(i,j)=cf*neffmx(i,j)*coef_ile ! cat betamx(i,j)=min(betamx(i,j),tob_ile) ! version Jorge ! ile : frott maxi=20, mini=2 ! betamx(i,j)=min(betamx(i,j),betailemax) ! betamx(i,j)=max(betamx(i,j),betailemin) else if (flotmx(i,j)) then ! flottant ou ile betamx(i,j)=0. else ! grounded, SIA betamx(i,j)=tostick ! frottement glace posee (1 bar) endif end do end do do j=1,ny ! le noeud 1 n'est pas attribue do i=1,nx if ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) gzmy(i,j)=.true. ! calcul du frottement basal (ce bloc etait avant dans neffect) if (cotemy(i,j)) then ! point cotier betamy(i,j)=cf*neffmy(i,j)*0.1 ! cat betamy(i,j)=min(betamy(i,j),tobmax) ! version Jorge ! point cotier: frott maxi=40, mini=4 ! betamy(i,j)=min(betamy(i,j),betacotemax) ! betamy(i,j)=max(betamy(i,j),betacotemin) ! version Jorge stream en zone sediment : else if ((fleuve_sedimy(i,j).and.gzmx(i,j).and.(.not.cotemx(i,j)))) then ! stream en zone sediment ! stream sediment: frott maxi=60, mini=6 !cdc betamy(i,j)=cf*neffmy(i,j)*3 betamy(i,j)=cf*neffmy(i,j)*coef_sedim ! betamy(i,j)=min(betamy(i,j),betasedmax) ! betamy(i,j)=max(betamy(i,j),betasedmin) else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then ! stream normal ! cat betamy(i,j)=tobmax ! cat betamy(i,j)=cf*neffmy(i,j)*20. ! cat betamy(i,j)=min(betamy(i,j),tobmax) ! cat betamy(i,j)=max(betamy(i,j),10.) ! version Jorge ! stream normal: frott maxi=80, mini=8 betamy(i,j)=cf*neffmy(i,j)*coef_gz ! betamy(i,j)=min(betamy(i,j),betastreamax) ! betamy(i,j)=max(betamy(i,j),betastreamin) ! cat if (cf*neffmy(i,j).gt.1500.) then ! cat gzmy(i,j)=.false. ! cat betamy(i,j)=tostick ! cat endif else if (ilemy(i,j)) then betamy(i,j)=cf*neffmy(i,j)*coef_ile ! cat betamy(i,j)=min(betamy(i,j),tob_ile) ! version Jorge ! ile : frott maxi=20, mini=2 ! betamy(i,j)=min(betamy(i,j),betailemax) ! betamy(i,j)=max(betamy(i,j),betailemin) else if (flotmy(i,j)) then ! flottant ou ile betamy(i,j)=0. else ! grounded, SIA betamy(i,j)=tostick ! frottement glace posee endif end do end do !print*,'dragging debug h_sedim(78,144),hwater,seuil_sedim,seuil_hwater',h_sedim(78,144),seuil_sedim,hwater(78,144),seuil_hwater,fleuve(78,144),gzmx(78,144) !print*,'dragging debug h_sedim(90,157),hwater,seuil_sedim,seuil_hwater',h_sedim(90,157),seuil_sedim,hwater(90,157),seuil_hwater,fleuve(90,157),gzmx(90,157),betamx(90,157),betamx(78,144) where (fleuve(:,:)) debug_3D(:,:,1)=1 elsewhere (flot(:,:)) debug_3D(:,:,1)=2 elsewhere debug_3D(:,:,1)=0 endwhere where (cote(:,:)) debug_3D(:,:,2)=1 elsewhere debug_3D(:,:,2)=0 endwhere where (fleuvemx(:,:)) debug_3D(:,:,3)=1 elsewhere debug_3D(:,:,3)=0 endwhere where (flotmx(:,:)) debug_3D(:,:,3)=-1 endwhere !_____________________ where (fleuvemy(:,:)) debug_3D(:,:,4)=1 elsewhere debug_3D(:,:,4)=0 endwhere where (flotmy(:,:)) debug_3D(:,:,4)=-1 end where debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) !print*,'point (79,145)',hwater(79,145),h_sedim(79,145),flot(79,145),fleuve_sedim(79,145),fleuve(79,145),betamx(79,145),fleuve_sedimx(79,145),gzmx(79,145) return end subroutine dragging end module dragging_hwat_sedim