!< \file beta_iter_vitbil_mod.f90 !< Module pour affiner le beta a partir d'une carte vitbil sur les noeuds centres. !< A appeler depuis steps_time_loop !< a utiliser avec beta_prescr !> \ beta_iter_vitbil !! \author ... !! \date ... !! @note Used module !! @note - use module3D_phy !< module beta_iter_vitbil_mod use module3d_phy use interface_input use netcdf use io_netcdf_grisli implicit none real*8, dimension(:,:), pointer :: tab => null() !< tableau 2d real ecrit dans le fichier pour lect ncdf 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 :: 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) real, dimension(nx,ny) :: Umag_vitbil !< vitesse de bilan real, dimension(nx,ny) :: Uslid_vitbil !< vitesse de bilan partie glissement real, dimension(nx,ny) :: Umag_direct !< vitesse moyenne centree calculee par GRISLI real, dimension(nx,ny) :: Uslmag_direct !< vitesse glissement centree calculee par GRISLI real, dimension(nx,ny) :: Udefmag_direct !< vitesse deformation centree calculee par GRISLI real, dimension(nx,ny) :: driving_stress !< driving stress character(len=100) :: Umag_bil_file !< fichier qui contient les donnees real :: time_iter !< temps de debut des iterations integer :: nb_iter_vitbil !< nombre d'iteratiosn avant de changer de pas de temps real :: coef_iter_vitbil !< coefficient <= 1 (rapport vitesses) real :: J_Umag !< estimateur ecart entre vitesses real :: J_Udef !< estimateur vitesse regions sans glissement integer :: nb_umag !< nombre de points calcul J_Umag integer :: nb_Udef !< nombre de points calcul J_Umag 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 ! call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'') write(6,*) beta_c_file call Read_Ncdf_var('z',beta_c_file,tab) drag_centre(:,:) = tab(:,:) ! 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(flot(:,:)) drag_centre(:,:) = 0. end where 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 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) ! 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 init_beta_iter_vitbil return end subroutine init_dragging !________________________________________________________________________________ !----------------------------------------------------------------------------------------- !< subroutine : init_beta_iter_vitbil !< Cette routine fait l'initialisation de la procedure beta_iter_vitbil subroutine init_beta_iter_vitbil namelist/beta_iter_vitbil/time_iter,nb_iter_vitbil,coef_iter_vitbil,Umag_bil_file rewind(num_param) ! pour revenir au debut du fichier param_list.dat 428 format(A) read(num_param,beta_iter_vitbil) write(num_rep_42,428) '!___________________________________________________________' write(num_rep_42,428) '! module beta_iter_vitbil' write(num_rep_42,beta_iter_vitbil) write(num_rep_42,428) '! ' write(num_rep_42,428) '! Umag_bil_file : fichier qui contient les donnees Umag_vitbil' write(num_rep_42,428) '! time_iter : temps de debut des iterations ' write(num_rep_42,428) '! nb_iter_vitbil: nombre diteratation a faire avant de changer de pas de temps ' write(num_rep_42,*) '! coef_iter_vitbil : coefficient pour le rapport vitesses ' write(num_rep_42,428) '!___________________________________________________________' time_iter = time_iter + tbegin ! time_iter est le temps de relaxation (de l'ordre de 5 ans) ! ajouter tbegin pour ne pas dependre du temps de depart. Umag_bil_file = trim(dirnameinp)//trim(Umag_bil_file) write(6,*) Umag_bil_file ! call lect_input(1,'Umag_vitbil',1,Umag_vitbil,Umag_bil_file,'') call Read_Ncdf_var('Umag_vitbil',Umag_bil_file,tab) Umag_vitbil(:,:) = tab(:,:) debug_3D(:,:,58) = Umag_vitbil(:,:) end subroutine init_beta_iter_vitbil !________________________________________________________________________________ !> 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') 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 !______________________________________________________________________________________ !>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_c2stag !! 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_iter_vitbil !< fait des iterations pour affiner la valeur de beta_centre pour s'ajuster sur vitbil !----------------------------------------------------------------------------------------- subroutine beta_iter_vitbil ! calcule la vitesse centree venant du calcul direct Umag_direct(:,:) = (((uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))**2+ & ! moyenne (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5 Uslmag_direct(:,:) = (((ux(:,:,nz)+eoshift(ux(:,:,nz),shift=1,boundary=0.,dim=1))**2+ & ! glissement (uy(:,:,nz)+eoshift(uy(:,:,nz),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5 Udefmag_direct(:,:) = max( (Umag_direct(:,:)- Uslmag_direct(:,:)) , 0.) ! deformation ! guess du beta dans les zones où GRISLI a une vitesse de glissement nulle call slope_surf driving_stress(:,:) = rog * H(:,:) * slope(:,:) ! driving stress Uslid_vitbil(:,:) = Umag_vitbil(:,:) - Udefmag_direct(:,:) ! calcule la vitesse de glissement "bilan" ! Pas de glissement dans GRISLI where ((Uslmag_direct(:,:).le.0.).and.(Uslid_vitbil(:,:).gt.1.e-3)) ! pas de glissement dans GRISLI beta_centre(:,:) = driving_stress(:,:) / Uslid_vitbil(:,:) ! mais du glissement "bilan" end where where ((Uslmag_direct(:,:).le.0.).and.(Uslid_vitbil(:,:).le.1.e-3)) ! ni l'un ni l'autre Uslid_vitbil(:,:) = 0. beta_centre(:,:) = beta_limgz end where ! Glissement dans GRISLI where (Uslid_vitbil(:,:).le.1.e-3) ! la vitesse de bilan est deja trop rapide Uslid_vitbil(:,:) = 0. beta_centre(:,:) = beta_limgz elsewhere (Uslmag_direct(:,:).gt.0.) beta_centre(:,:) = beta_centre(:,:) * coef_iter_vitbil * Uslmag_direct(:,:) / Uslid_vitbil(:,:) end where where(flot(:,:)) beta_centre(:,:) = 0. end where beta_centre(:,:) = min(beta_centre(:,:),beta_limgz) call beta_c2stag(nx,ny,betamx,betamy,beta_centre) ! redistribue sur les mailles alternees call beta_min_value(beta_min) ! valeur mini que peut avoir beta (en Pa an m-1) debug_3D(:,:,59) = betamx(:,:) debug_3D(:,:,60) = beta_centre(:,:) ! estimateurs J_Umag = 0. J_Udef = 0. nb_umag = 0 nb_udef = 0 do j=2,ny-1 do i=2,nx-1 if (.not.flot(i,j)) then ! calcul sur les points poses nb_umag = nb_umag + 1 J_Umag = (Umag_direct(i,j) - Umag_vitbil(i,j))**2. + J_Umag if (Uslid_vitbil(i,j).le.1.e-3) then ! calcul sur la partie deformation suelement nb_udef = nb_udef + 1 J_Udef = (Umag_direct(i,j) - Umag_vitbil(i,j))**2. + J_Udef end if end if end do end do J_Umag = (J_Umag**0.5) / nb_umag J_Udef = (J_Udef**0.5) / nb_udef drag_centre(:,:) = beta_centre(:,:) drag_mx(:,:) = betamx(:,:) drag_my(:,:) = betamy(:,:) i=319 j=113 write(266,'(2(i0,1x),3x,9(es12.4,1x))') nb_umag,nb_udef,time,J_Umag,J_Udef,Umag_direct(i,j),Umag_vitbil(i,j),Uslmag_direct(i,j),Uslid_vitbil(i,j),beta_centre(i,j),betamx(i,j),drag_centre(i,j) end subroutine beta_iter_vitbil !----------------------------------------------------------------------------------------- end module beta_iter_vitbil_mod