!< \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, 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 real :: time_iter_end !< temps de fin des iterations real :: time_reiter !< nbr annees entre 2 iterations calcul beta 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_beta_iter_vitbil !< Cette routine fait l'initialisation de la procedure beta_iter_vitbil subroutine init_beta_iter_vitbil implicit none real*8, dimension(:,:), pointer :: tab => null() !< tableau 2d real ecrit dans le fichier pour lect ncdf namelist/beta_iter_vitbil/time_iter,nb_iter_vitbil,coef_iter_vitbil,Umag_bil_file,time_reiter 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,*) '! time_reiter : nombre d annees entre 2 iterations' 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. time_iter_end = time_iter + 1. ! 20. ans en version std Cat et Seb print*,'debug init_beta_iter_vitbil time_iter=', time_iter, time_iter_end, time 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: 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(m_iter) implicit none integer :: m_iter ! indice bloucle iteration ! calcul des vitesses cibles : !if ((time.eq.time_iter).and.(time.gt.time_reiter).and.(m_iter.eq.1)) then if ((abs(time-time_iter).lt.dtmin).and.(time.gt.time_reiter).and.(m_iter.eq.1)) then 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 Umag_vitbil(:,:)= min(Umag_direct(:,:) * min(1.e3,H(:,:)/H0(:,:)), 5000.) print*,'*****debug beta_iter_vitbil calcul de Umag_vitbil',time,m_iter do j=1,ny do i=1,nx write(266,*) Umag_vitbil(i,j) enddo enddo endif ! 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(:,:) end subroutine beta_iter_vitbil !----------------------------------------------------------------------------------------- end module beta_iter_vitbil_mod