!> \file dragging_prescr_beta_mod.f90 !! Module pour calculer le beta a partir d'une carte sur les noeuds centres. !< !> \namespace dragging_prescr_beta !! Calcule le beta a partir de vitesses de bilan !! @note Il faut partir d'un cptr !! \author ... !! \date ... !! @note Used module !! @note - use module3D_phy !< module dragging_prescr_beta use module3d_phy use interface_input implicit none real, dimension(nx,ny) :: Vcol_x !< uniquement pour compatibilite real, dimension(nx,ny) :: Vcol_y !< real, dimension(nx,ny) :: Vsl_x !< real, dimension(nx,ny) :: Vsl_y !< real, dimension(nx,ny) :: drag_centre !< beta on major node (average) real, dimension(nx,ny) :: betamoy_x !< pour faire les moyennes glissantes real, dimension(nx,ny) :: betamoy_y !< integer, dimension(nx,ny) :: nbvoisx !< integer, dimension(nx,ny) :: nbvoisy !< real :: beta_limgz !< when beta gt beta_limgz -> not gzmx real :: beta_min !< minimum value of beta real :: beta_mult !< coefficient of beta field integer :: ill,jll,nmoy real :: maxi !< calcul correction dS a la grounding line real :: mini logical :: corr_def = .False. !< for deformation correction (compatibility) contains !----------------------------------------------------------------------------------------- !> SUBROUTINE: init_dragging !! Cette routine fait l'initialisation du dragging. !< subroutine init_dragging ! Cette routine fait l'initialisation du dragging. ! nouvelle version : lit les fichiers nc implicit none character(len=100) :: beta_c_file ! beta on centered grid character(len=100) :: betamx_file ! beta on staggered grid mx character(len=100) :: betamy_file ! beta on staggered grid mx namelist/beta_prescr/beta_c_file,beta_limgz,beta_min,beta_mult if (itracebug.eq.1) call tracebug(' Prescr_beta subroutine init_dragging') iter_beta=0 ! lecture des parametres du run !-------------------------------------------------------------------- rewind(num_param) ! pour revenir au debut du fichier param_list.dat 428 format(A) read(num_param,beta_prescr) write(num_rep_42,428) '!___________________________________________________________' write(num_rep_42,428) '! read beta on centered grid' write(num_rep_42,beta_prescr) write(num_rep_42,428) '! beta_file : nom des fichier qui contiennent les beta_c' write(num_rep_42,428) '! above beta_limgz, gzmx is false' write(num_rep_42,428) '! if grounded, beta_min minimum value ' write(num_rep_42,428) '! beta_mult : coefficient just after reading' write(num_rep_42,*) write(num_rep_42,428) '!___________________________________________________________' beta_c_file = trim(dirnameinp)//trim(beta_c_file) betamx_file = trim(dirnameinp)//trim(betamx_file) betamy_file = trim(dirnameinp)//trim(betamy_file) betamax = beta_limgz ! read the beta value on centered and staggered grids ! call lect_input(1,'beta_c',1,drag_centre,beta_c_file,trim(dirnameinp)//trim(runname)//'.nc') ! call lect_input(1,'betamx',1,drag_mx,betamx_file,trim(dirnameinp)//trim(runname)//'.nc') ! call lect_input(1,'betamy',1,drag_my,betamy_file,trim(dirnameinp)//trim(runname)//'.nc') call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'') ! multiplication par le coefficient beta_mult drag_centre(:,:) = drag_centre(:,:)*beta_mult where (mk_init(:,:).eq.-2) ! iles a probleme drag_centre(:,:) = 2.* abs(beta_limgz) end where where ((H(:,:).lt.0.1).and.(Bsoc(:,:).gt.0)) ! zones actuellement non couvertes drag_centre(:,:) = 1000. ! valeur relativement faible qui evite end where ! les freinages ! call lect_input(1,'betamx',1,drag_mx,betamx_file,'') ! call lect_input(1,'betamy',1,drag_my,betamy_file,'') call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre) ! redistribute on staggered grid if (beta_limgz.gt.0.) then beta_centre(:,:) = drag_centre(:,:) betamx(:,:) = drag_mx(:,:) betamy(:,:) = drag_my(:,:) else if (beta_limgz.lt.0.) then ! attention sans doute pour plastic beta_centre(:,:) = - beta_limgz betamx(:,:) = - beta_limgz betamy(:,:) = - beta_limgz drag_centre(:,:) = - beta_limgz drag_mx(:,:) = - beta_limgz drag_my(:,:) = - beta_limgz beta_limgz = 1. end if ! divers traitements qui peuvent etre appliques a beta : mis sous forme de routine ! call beta_plateau_slope ! give value beta_limgz to places where slope is >1.e-2 ! call beta_stag2c ! average the staggered values on central ones ! call beta_plateau_S ! give value beta_limgz to places where S > 2500 m call beta_min_value(beta_min) ! valeur mini que peut avoir beta (en Pa an m-1) ! on peut envisager une valeur ~ 10 ! rappel : beta doit être positif. ! call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre) ! redistribute on staggered grid ! mstream attribution !________________________ ! mstream -> regions where streaming is allowed ! drag_mx -> value of beta if an ice stream is allowed ! betamx -> value of dragging (may depend on other conditions) ! to ensure having some "Dirichlet" points, there must be at least ! a mstream=0 in every island ! un tour de lissage ! call beta_stag2c(nx,ny,drag_mx,drag_my,drag_centre) ! average staggered ! ! on central ! call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre) ! centered distributed ! ! on staggered ! beta_centre(:,:) = drag_centre(:,:) ! surtout pour sorties where (drag_mx(:,:).ge.beta_limgz) mstream_mx(:,:) = 0 betamx(:,:) = beta_limgz drag_mx(:,:) = beta_limgz elsewhere mstream_mx(:,:) = 1 betamx(:,:) = drag_mx(:,:) endwhere where (drag_my(:,:).ge.beta_limgz) mstream_my(:,:) = 0 betamy(:,:) = beta_limgz drag_my(:,:) = beta_limgz elsewhere mstream_my(:,:) = 1 betamy(:,:) = drag_my(:,:) endwhere if (itracebug.eq.1) call tracebug(' Prescr_beta apres lecture') mstream_mx(:,:) = 0 mstream_my(:,:) = 0 call gzm_beta_prescr ! call write_datfile2(nx,ny,betamx,betamy,'beta-LBq15-08.dat') ! call make_gaussienne(n_g,sigma_beta,M_gauss) return end subroutine init_dragging !________________________________________________________________________________ !> SUBROUTINE: dragging !! Defini la localisation des streams et le frottement basal !< !------------------------------------------------------------------------- subroutine dragging ! defini la localisation des streams et le frottement basal implicit none if (itracebug.eq.1) call tracebug(' beta_prescr subroutine dragging') ! tentative d'iteration ! call beta_stag2c(nx,ny,drag_mx,drag_my,drag_centre) ! average staggered ! on central ! call iter_gauss(n_g,nx,ny,M_gauss,-1.,drag_centre,corrstream_centre,hdot) ! correct central ! call beta_c2stag(nx,ny,drag_mx,drag_my,corrstream_centre) ! centered distributed ! on staggered ! drag_centre(:,:) = corrstream_centre(:,:) ! beta_centre(:,:) = drag_centre(:,:) ! surtout pour sorties ! pour l'instant pas d'autre condition que mstream where (mstream_mx(:,:).eq.1) betamx(:,:) = drag_mx(:,:) elsewhere betamx(:,:)= beta_limgz end where where (mstream_my(:,:).eq.1) betamy(:,:) = drag_my(:,:) elsewhere betamy(:,:)= beta_limgz end where betamx(:,:)=drag_mx(:,:) betamy(:,:)=drag_my(:,:) call gzm_beta_prescr ! determine gz et flgz et met a 0 le beta des points flottants return end subroutine dragging !---------------------------------------------------------------------------------- !>SUBROUTINE: gzm_beta_prescr !! Calcul de gzmx !< subroutine gzm_beta_prescr ! attribue flgzmx et gzmx (et my) logical :: test_Tx ! test if there is a basal melting point in the surrounding logical :: test_Ty ! test if there is a basal melting point in the surrounding logical :: test_beta_x ! test if there is a low drag point in the surrounding logical :: test_beta_y ! test if there is a low drag point in the surrounding real :: beta_forc_fleuve ! below this value -> ice stream ! beta_forc_fleuve = 5.e3 beta_forc_fleuve = beta_limgz ! calcul de gzmx gzmx(:,:)=.true. gzmy(:,:)=.true. flgzmx(:,:)=.false. flgzmy(:,:)=.false. ! points flottants : flgzmx mais pas gzmx do j=2,ny do i=2,nx if (flot(i,j).and.(flot(i-1,j))) then flotmx(i,j)=.true. flgzmx(i,j)=.true. gzmx(i,j)=.false. betamx(i,j)=0. end if if (flot(i,j).and.(flot(i,j-1))) then flotmy(i,j)=.true. flgzmy(i,j)=.true. gzmy(i,j)=.false. betamy(i,j)=0. end if end do end do if (itracebug.eq.1) call tracebug(' apres criteres flot') if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330) !--------- autres criteres ! points poses gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz) ! Pas de calcul pour les points gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz) ! au fort beta ! critere (lache) sur la temperature do j = 2, ny-1 do i =2, nx-1 ! test s'il y a un point tempere dans les environs test_Tx = any( ibase( i-1:i , j-1:j+1 ).gt.1) test_Ty = any( ibase( i-1:i+1 , j-1:j ).gt.1) ! test s'il y a un point low drag dans les environs test_beta_x = any( drag_centre( i-1:i , j-1:j+1 ) .lt. beta_forc_fleuve ) test_beta_y = any( drag_centre( i-1:i+1 , j-1:j ) .lt. beta_forc_fleuve ) ! critere pour gzmx ! Attention : les deux lignes suivants sont en commentaire pour voir l'effet ! gzmx(i,j) = gzmx(i,j) .and. (test_Tx.or.test_beta_x) ! gzmy(i,j) = gzmy(i,j) .and. (test_Ty.or.test_beta_y) end do end do if (itracebug.eq.1) call tracebug(' apres autres criteres ') if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330) flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:) flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:) where ((.not.flotmx(:,:)).and.(.not.gzmx(:,:))) betamx(:,:) = beta_limgz end where where ((.not.flotmy(:,:)).and.(.not.gzmy(:,:))) betamy(:,:) = beta_limgz end where end subroutine gzm_beta_prescr !---------------------------------------------------------------------------------- ! Routines qui permettent des manip simples sur beta. !______________________________________________________________________________________ !>SUBROUTINE: beta_plateau_slope !! give value beta_limgz to places where slope is >1.e-2 !< subroutine beta_plateau_slope ! give value beta_limgz to places where slope is >1.e-2 call slope_surf where (abs(sdx(:,:)).gt.1.e-2) drag_mx(:,:) = beta_limgz end where where (abs(sdy(:,:)).gt.1.e-2) drag_my(:,:) = beta_limgz end where end subroutine beta_plateau_slope !______________________________________________________________________________________ !> SUBROUTINE beta_plateau_S !! give value beta_limgz to places where S > 2500 m !< subroutine beta_plateau_S ! give value beta_limgz to places where S > 2500 m where (S(:,:).gt.2500.) ! moyen simple de nettoyer l'interieur. WAIS pas touche drag_centre(:,:)=beta_limgz drag_mx(:,:)=beta_limgz drag_my(:,:)=beta_limgz end where end subroutine beta_plateau_S !______________________________________________________________________________________ !>SUBROUTINE: beta_min_value !! valeur mini que peut avoir beta (en Pa an m-1) !< subroutine beta_min_value(val_betamin) real :: val_betamin ! valeur mini que peut avoir beta (en Pa an m-1) drag_centre(:,:) = max(drag_centre(:,:),val_betamin) drag_mx(:,:) = max(drag_mx(:,:),val_betamin) drag_my(:,:) = max(drag_my(:,:),val_betamin) where(flot(:,:)) drag_centre(:,:) = 0. end where end subroutine beta_min_value !______________________________________________________________________________________ !>SUBROUTINE: beta_stag2c !! average the staggered values on central ones !< subroutine beta_stag2c(nxx,nyy,R_mx,R_my,R_centre) ! average the staggered values on central ones implicit none integer :: nxx,nyy ! dimension des tableaux real, dimension(nxx,nyy), intent(in) :: R_mx ! valeur du beta sur maille mx real, dimension(nxx,nyy), intent(in) :: R_my ! valeur du beta sur maille my real, dimension(nxx,nyy), intent(out) :: R_centre ! valeur du beta sur maille centree R_centre(:,:)=0. do j=2,ny-1 do i=2,nx-1 R_centre(i,j)= ((R_mx(i,j)+R_mx(i+1,j)) & + (R_my(i,j)+R_my(i,j+1)))*0.25 end do end do end subroutine beta_stag2c !______________________________________________________________________________________ !>SUBROUTINE: beta_c2stag_min !! redistribute central values on staggered grid !< subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid implicit none integer :: nxx,nyy ! dimension des tableaux real, dimension(nxx,nyy), intent(out) :: R_mx ! valeur du beta sur maille mx real, dimension(nxx,nyy), intent(out) :: R_my ! valeur du beta sur maille my real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree do j=2,nyy do i=2,nxx R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5 R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5 end do end do end subroutine beta_c2stag !______________________________________________________________________________________ !>SUBROUTINE: beta_c2stag_min !! redistribute central values on staggered grid !< subroutine beta_c2stag_min(nxx,nyy,R_mx,R_my,R_centre) ! redistribute central values on staggered grid implicit none integer :: nxx,nyy ! dimension des tableaux real, dimension(nxx,nyy), intent(inout):: R_mx ! valeur du beta sur maille mx real, dimension(nxx,nyy), intent(inout):: R_my ! valeur du beta sur maille my real, dimension(nxx,nyy), intent(in) :: R_centre ! valeur du beta sur maille centree do j=2,nyy do i=2,nxx R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5) R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5) end do end do end subroutine beta_c2stag_min !______________________________________________________________________________________ !>SUBROUTINE: make_gaussienne !! cree une matrice gaussienne normalisee a 1 !< subroutine make_gaussienne(n,sigma,M) implicit none integer,intent(in) :: n ! demi taille du tableau gaussienne real, intent(in) :: sigma ! sigma gaussienne exprime en dx real,dimension(-n:n,-n:n),intent(inout) :: M ! tableau gaussienne ! variables locales real :: dist2 ! distance au centre (au carre) real :: som ! pour calibration real :: sigma2 ! pour gaussienne som = 0. sigma2 = sigma * sigma * 2 do j = -n, n do i = -n, n dist2 = (i*i+j*j) M(i,j) = exp(-dist2/sigma2) som=som+M(i,j) end do end do M(:,:) = M(:,:) / som ! do j = -n, n ! write(6,'(5(f0.4,1x))') (M(i,j),i=-n,n) ! end do end subroutine make_gaussienne !______________________________________________________________________________________ !>SUBROUTINE: iter_gauss !! utilise hdot pour faire les iterations sur beta_centre !< subroutine iter_gauss(n,nxx,nyy,M,coef,beta_in,beta_out,Ecart) integer, intent(in) :: n ! taille matrice gaussienne integer, intent(in) :: nxx,nyy ! taille du domaine real,dimension(-n:n,-n:n),intent(in) :: M ! matrice gaussienne real, intent(in) :: coef ! coefficient multiplicateur real,dimension(nx,ny), intent(in) :: beta_in ! la valeur initiale real,dimension(nx,ny), intent(out) :: beta_out ! la valeur apres iteration real,dimension(nx,ny), intent(in) :: Ecart ! l'ecart a l'obs ! variables locales integer :: ii,jj do j= 1 + n, nyy - n do i = 1 + n, nxx - n beta_out(i,j) = beta_in(i,j) do jj = -n, n do ii = -n, n beta_out(i,j) = beta_out(i,j) + coef * M(ii,jj) * ecart(i+ii,j+jj) end do end do end do end do beta_out(:,:) = max (beta_out(:,:),0.) end subroutine iter_gauss end module dragging_prescr_beta !!$! add the nodes of topographical instability (on major nodes) ! a garder eventuellement pour runs instabilite !!$beta_file = trim(dirnameinp)//'sensib_H_Gael-15km.dat' !!$call lect_datfile(nx,ny,drag_instab_topo,1,beta_file) !!$ !!$where (drag_instab_topo(:,:).gt.0.) !!$ drag_centre(:,:)=min(drag_centre(:,:),400.) !!$endwhere