!> \file dragging_neff_contmaj_mod.f90 !! Defini les zones de stream avec 3 criteres !< !> \namespace dragging_neff_contmaj !! Defini les zones de stream avec 3 criteres !! \author ... !! \date ... !! @note Les trois criteres sont: !! @note * un critere sur la pression effective !! @note * un critere de continuite depuis la cote !! @note * un critere sur la courbure du socle (si negatif, vallees) !! @note Used module !! @note - use module3D_phy !! @note - use param_phy_mod !< module dragging_neff_slope ! Defini les zones de stream avec : ! * un critere sur la pression effective ! * un critere de continuite depuis la cote ! * un critere sur la courbure du socle (si negatif, vallees) use module3d_phy use param_phy_mod use interface_input use io_netcdf_grisli use fake_beta_iter_vitbil_mod implicit none logical,dimension(nx,ny) :: fleuve logical,dimension(nx,ny) :: cote logical,dimension(nx,ny) :: slowssamx ! not actual stream, but ssa used as Ub logical,dimension(nx,ny) :: slowssamy ! not actual stream, but ssa used as Ub logical,dimension(nx,ny) :: slowssa ! not actual stream, but ssa used as Ub real,dimension(nx,ny) :: hires_slope ! slope comupted on a high resolution topography character(len=100) :: slope_fich ! fichier grille character(len=100) :: file_ncdf !< fichier netcdf issue des fichiers .dat 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 :: betamin ! betamin : (Pa) frottement mini sous les streams real :: tob_ile ! pour les iles real :: cry_lim=50. ! courbure limite pour le suivi des fleuves real,dimension(nx,ny) :: neff ! pression effective noeuds majeurs real :: seuil_neff ! seuil sur la pression effective pour avoir glissement en zone sediment real :: coef_gz ! coef frottement zones stream std (10) real :: coef_ile ! coef frottement zones iles (0.1) 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 contains !------------------------------------------------------------------------------- !> SUBROUTINE: init_dragging !! Cette routine fait l'initialisation du dragging. !> subroutine init_dragging ! Cette routine fait l'initialisation du dragging. implicit none ! local variables, defining the rugosity-enhanced flow real :: expo_slope real :: pente_min, pente_max namelist/drag_neff_slope/cf,betamax,betamin,toblim,tostick,seuil_neff,coef_gz,coef_ile,slope_fich,expo_slope,pente_min,pente_max if (itracebug.eq.1) call tracebug(' dragging avec neff et slope') ! formats pour les ecritures dans 42 428 format(A) ! lecture des parametres du run block drag neff !-------------------------------------------------------------------- rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,drag_neff_slope) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '&drag_neff_slope ! nom du bloc dragging neff slope' write(num_rep_42,*) write(num_rep_42,*) 'cf = ',cf write(num_rep_42,*) 'betamax = ', betamax write(num_rep_42,*) 'betamin = ', betamin write(num_rep_42,*) 'toblim = ', toblim write(num_rep_42,*) 'tostick = ', tostick write(num_rep_42,*) 'seuil_neff = ', seuil_neff write(num_rep_42,*) 'coef_gz = ', coef_gz write(num_rep_42,*) 'coef_ile = ', coef_ile write(num_rep_42,'(A,A)') 'slope_fich = ', slope_fich write(num_rep_42,*) 'expo_slope = ', expo_slope write(num_rep_42,*) 'pente_min = ', pente_min write(num_rep_42,*) 'pente_max = ', pente_max write(num_rep_42,*)'/' write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff' write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams ' write(num_rep_42,428) '! betamin : (Pa) frottement mini sous les streams ' write(num_rep_42,428) '! toblim : (Pa) pour les iles ' write(num_rep_42,428) '! tostick : (Pa) pour les points non flgzmx ' write(num_rep_42,428) '! seuil_neff (Pa) seuil sur la pression effective pour avoir glissement' 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,428) '! slope_fich : fichier pente bedrock' write(num_rep_42,428) '! expo_slope : exposant fonction pente' write(num_rep_42,428) '! pente_min : no flow enhancement below this' write(num_rep_42,428) '! pente_max : maximum flow enhancement above this' write(num_rep_42,*) inv_beta=0 tob_ile=betamax/2. moteurmax=toblim ! betamax_2d depends on sub-grid slopes. slope_fich=trim(dirnameinp)//trim(slope_fich) call lect_input(1,'slope',1,hires_slope(:,:),slope_fich,file_ncdf) ! from slopes, we create an index between 0 & 1 ! 1 is very mountainous, 0 is flat hires_slope(:,:)=(max(min(hires_slope(:,:),pente_max),pente_min)-pente_min)/(pente_max-pente_min) ! now we compute the actual betamax used by the remplimat routines ! /!\ the slope is used to modify the beta where we have a temperate base, ! but NO ice stream... -> Slow SSA zone (SSA used to compute Ub) betamax_2d(:,:)=max ( tostick * (1. - (1 - betamax / tostick ) * hires_slope(:,:)**(1./expo_slope)) , betamax ) !do j=1,ny ! do i=1,nx ! write(18745,*) betamax_2d (i,j) ! enddo !enddo !------------------------------------------------------------------- ! 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 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 !$ USE OMP_LIB ! les iles n'ont pas de condition neff mais ont la condition toblim ! (ce bloc etait avant dans flottab) !$OMP PARALLEL !$OMP DO 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 !$OMP 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. !$OMP WORKSHARE fleuvemx(:,:)=.false. fleuvemy(:,:)=.false. fleuve(:,:)=.false. cote(:,:)=.false. !$OMP END WORKSHARE ! detection des cotes !$OMP DO 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 !$OMP END DO ! calcul de neff (pression effective sur noeuds majeurs) !$OMP DO do j=1,ny-1 do i=1,nx-1 neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4 enddo enddo !$OMP END DO !aurel, for the borders: !$OMP WORKSHARE neff(:,ny)=1.e5 neff(nx,:)=1.e5 ! aurel, we add the neff threshold: where ((neff(:,:).le.seuil_neff).and.(.not.flot(:,:)).and.(H(:,:).gt.1.)) fleuve(:,:)=.true. !$OMP END WORKSHARE !$OMP 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 end do end do !$OMP END DO ! we look for the warm based points that will not be treated as stream (ub from SSA): !$OMP WORKSHARE slowssa(:,:)=.false. slowssamx(:,:)=.false. slowssamy(:,:)=.false. !$OMP END WORKSHARE !$OMP DO do j=1,ny do i=1,nx !if ((not(flot(i,j))).and.(hwater(i,j).gt.0.1)) slowssa(i,j)=.true. if ((.not.(flot(i,j))).and.(ibase(i,j).ne.1).and.(H(i,j).gt.1.)) slowssa(i,j)=.true. end do end do !$OMP END DO !$OMP DO do j=1,ny-1 do i=1,nx-1 if (slowssa(i,j)) then slowssamx(i,j)=.true. slowssamx(i+1,j)=.true. slowssamy(i,j)=.true. slowssamy(i,j+1)=.true. end if end do end do !$OMP END DO !$OMP DO do j=1,ny do i=1,nx ! the actual streams and the warm based points are gzmx now: if ( ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) .or. ((.not.ilemx(i,j)).and.(slowssamx(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) betamx(i,j)=min(betamx(i,j),betamax) else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then ! tous les points SSA if (fleuvemx(i,j)) then ! the actual streams betamx(i,j)=cf*neffmx(i,j)*coef_gz betamx(i,j)=min(betamx(i,j),betamax) betamx(i,j)=max(betamx(i,j),betamin) if (cf*neffmx(i,j).gt.betamax*2.) then ! a stream that's becoming grounded... if (slowssamx(i,j)) then betamx(i,j)=betamax_2d(i,j) else gzmx(i,j)=.false. betamx(i,j)=tostick endif endif else ! not an actual stream, SSA is used to compute Ub betamx(i,j)=betamax_2d(i,j) endif else if (ilemx(i,j)) then betamx(i,j)=cf*neffmx(i,j)*coef_ile betamx(i,j)=min(betamx(i,j),tob_ile) 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 !$OMP END DO !$OMP DO do j=1,ny do i=1,nx ! the actual streams and the warm based points are gzmx now: if ( ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) .or. ((.not.ilemy(i,j)).and.(slowssamy(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) betamy(i,j)=min(betamy(i,j),betamax) else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then ! tous les points SSA if (fleuvemy(i,j)) then ! the actual streams betamy(i,j)=cf*neffmy(i,j)*coef_gz betamy(i,j)=min(betamy(i,j),betamax) betamy(i,j)=max(betamy(i,j),betamin) if (cf*neffmy(i,j).gt.betamax*2.) then ! a stream that's becoming grounded... if (slowssamy(i,j)) then betamy(i,j)=betamax_2d(i,j) else gzmy(i,j)=.false. betamy(i,j)=tostick endif endif else ! not an actual stream, SSA is used to compute Ub betamy(i,j)=betamax_2d(i,j) endif else if (ilemy(i,j)) then betamy(i,j)=cf*neffmy(i,j)*coef_ile betamy(i,j)=min(betamy(i,j),tob_ile) else if (flotmy(i,j)) then ! flottant ou ile betamy(i,j)=0. else ! grounded, SIA betamy(i,j)=tostick ! frottement glace posee (1 bar) endif end do end do !$OMP END DO !$OMP DO do j=2,ny-1 do i=2,nx-1 beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) & + (betamy(i,j)+betamy(i,j+1)))*0.25 end do end do !$OMP END DO !$OMP END PARALLEL !~ 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) return end subroutine dragging end module dragging_neff_slope