!> \file spinup_mod.f90 !! This module allows various spinup tasks mostly based on balance velocities. !< !> \namespace spinup_vitbil !! This module allows various spinup tasks mostly based on balance velocities. !! @note Could be transformed to use observed surface velocities !! @note must be used as a dragging module !! @note in module "choix" other dragging module must be switched off !! \author ... !! \date ... !! @note Used modules: !! @note - use module3D_phy !! @note - use param_phy_mod !! @note - use deformation_mod_2lois !< module spinup_vitbil use module3D_phy use param_phy_mod use deformation_mod_2lois use interface_input implicit none real, dimension(nx,ny) :: Vcol_x !< vertically averaged velocity x direction (balance) real, dimension(nx,ny) :: Vcol_y !< vertically averaged velocity y direction (balance) real, dimension(nx,ny) :: Ux_deformation !< vertically averaged deformation velocity x direction real, dimension(nx,ny) :: Uy_deformation !< vertically averaged deformation velocity y direction real, dimension(nx,ny) :: coef_defmx !< rescaling coefficient of Sa_mx and s2a_mx real, dimension(nx,ny) :: coef_defmy !< rescaling coefficient of Sa_my and s2a_my real, dimension(nx,ny) :: init_coefmx !< rescalling coefficient before limitation real, dimension(nx,ny) :: init_coefmy !< rescalling coefficient before limitation ! compatibilites avec remplimat logical :: corr_def = .False. !< for deformation correction (compatibility) real, dimension(nx,ny) :: Vsl_x !< real, dimension(nx,ny) :: Vsl_y !< contains !-------------------------------------------------------------------------------- !> SUBROUTINE: init_spinup !! Initialize spinup !< subroutine init_spinup namelist/spinup/ispinup ! put here which type of spinup ! ispinup = 0 ! run standard ! ispinup = 1 ! ispinup = 1 -> saute icethick3, diagnoshelf, ! diffusiv, isostasie pour equilibre temperature ! en prenant les vitesses calculees aux premiers ! pas de temps ! ispinup = 2 ! ispinup = 2 -> fait la conservation ! ! de la masse avec les vitesses de bilan ! ispinup = 3 ! fait le calcul des temperatures ! ! avec les vitesses de bilan ! lecture des parametres rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,spinup) write(num_rep_42,*)'!___________________________________________________________' write(num_rep_42,*)'! spinup module spinup_vitbil ' write(num_rep_42,spinup) ! pour ecrire les parametres lus write(num_rep_42,*) write(num_rep_42,*) write(num_rep_42,*)'! ispinup = 0 run standard voir no_spinup' write(num_rep_42,*)'! ispinup = 1 equilibre temperature avec vitesses grisli voir no_spinup' write(num_rep_42,*) write(num_rep_42,*)'! ispinup >1 use spinup_vitbil' write(num_rep_42,*)'! ispinup = 2 conservation de la masse avec vitesses bilan ' write(num_rep_42,*)'! ispinup = 3 equilibre temperature avec vitesses bilan' write(num_rep_42,*) if (ispinup.eq.0) then write(6,*)' ispinup = 0 et ispinup = 1 doivent etre appeles par no_spinup ' write(6,*) 'et il faut rajouter un dragging' stop endif call lect_vitbil if (itracebug.eq.1) call tracebug(' fin routine init_spinup de spinup_vitbil') end subroutine init_spinup subroutine lect_vitbil character(len=100) :: balance_Ux_file ! balance velocity file Ux character(len=100) :: balance_Uy_file ! balance velocity file Uy namelist/vitbil_upwind/balance_Ux_file, balance_Uy_file if (itracebug.eq.1) call tracebug(' Subroutine lect_vitbil') ! lecture des parametres du run !-------------------------------------------------------------------- rewind(num_param) ! pour revenir au debut du fichier param_list.dat 428 format(A) read(num_param,vitbil_upwind) write(num_rep_42,428) '!___________________________________________________________' write(num_rep_42,428) '!read balance velocities on staggered grid' write(num_rep_42,vitbil_upwind) write(num_rep_42,428) '! balance_Ux_file : nom du fichier qui contient les vit. bilan Ux' write(num_rep_42,428) '! balance_Uy_file : nom du fichier qui contient les vit. bilan Uy' write(num_rep_42,*) ! read balance velocities on staggered nodes balance_Ux_file = trim(dirnameinp)//trim(balance_Ux_file) balance_Uy_file = trim(dirnameinp)//trim(balance_Uy_file) call lect_input(1,'Vcol_x',1,Vcol_x,balance_Ux_file,"") call lect_input(1,'Vcol_y',1,Vcol_y,balance_Uy_file,"") ! call lect_input(2,'Vcol',1,Vcol_x,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc') !read Vcol_x ! call lect_input(2,'Vcol',2,Vcol_y,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc') !read Vcol_y !call lect_datfile(nx,ny,Vcol_x,1,balance_vel_file) ! read Vcol_x !call lect_datfile(nx,ny,Vcol_y,2,balance_vel_file) ! read Vcol_y debug_3D(:,:,59)=Vcol_x(:,:) debug_3D(:,:,60)=Vcol_y(:,:) end subroutine lect_vitbil !____________________________________________________________________________________ !> SUBROUTINE: init_dragging !! @ Cette routine fait l'initialisation du dragging. !< subroutine init_dragging ! Cette routine fait l'initialisation du dragging. implicit none mstream_mx(:,:)=0 ! pas de dragging a l'initialisation mstream_my(:,:)=0 ! coefficient permettant de modifier le basal drag. drag_mx(:,:)=1. drag_my(:,:)=1. gzmx(:,:)=.false. gzmy(:,:)=.false. flgzmx(:,:)=flotmx(:,:) flgzmy(:,:)=flotmy(:,:) end subroutine init_dragging subroutine dragging end subroutine dragging !---------------------------------------------------------------------- subroutine force_balance_vel if (itracebug.eq.1) call tracebug(' Subroutine force_balance_vel') Uxbar(:,:)=Vcol_x(:,:) Uybar(:,:)=Vcol_y(:,:) debug_3D(:,:,59)=Uxbar(:,:) debug_3D(:,:,60)=Uybar(:,:) end subroutine force_balance_vel !> SUBROUTINE: calc_coef_vitbil !! @ Calibrate Sa and S2a to force velocity to be balance_velocity in a consistent way. !! @note Must be called after flowlaw and before velocities !< subroutine calc_coef_vitbil ddbx(:,:)=0. ddby(:,:)=0. if (itracebug.eq.1) call tracebug(' Subroutine calc_coef_vitbil') call slope_surf where (abs(sdx(:,:)).lt.1.e-6) sdx(:,:)=-sign(1.e-6,Vcol_x(:,:)) end where where (abs(sdy(:,:)).lt.1.e-6) sdy(:,:)=-sign(1.e-6,Vcol_y(:,:)) end where call calc_ubar_def if (itracebug.eq.1) call tracebug(' apres calc_ubar_def') !---------------compute coef_defmx ----------------------------------- do j=1,ny do i=1,nx flotx: if (flotmx(i,j)) then coef_defmx(i,j) = 1. Uxbar(i,j) = Vcol_x(i,j) ! dmr below is a cast from a real*4 to logical*4 ! dmr cannot be implicit in gfortran ! dmr flgzmx(i,j) = Vcol_x(i,j) flgzmx(i,j) = (nint(Vcol_x(i,j)).ne.0) uxdef(i,j) = 0. Ubx(i,j) = Vcol_x(i,j) else coldx: if ((coefmxbmelt(i,j).le.0.).and.(.not.gzmx(i,j))) then !base froide !meme test que dans sliding Bindshadler ! sans doute trop fort ddbx(i,j) = 0. Ubx(i,j) = 0. if (abs(Ux_deformation(i,j)).gt.0.01) then ! calcule par calc_ubar_def coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j) else coef_defmx(i,j)=1. end if else ! base temperee if (abs(Ux_deformation(i,j)).gt.abs(Vcol_x(i,j))) then ! only deformation coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j) ddbx(i,j) = 0 Ubx(i,j) = 0. else coef_defmx(i,j)=1. ! no rescaling of deformation ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j) ! pb si directions opposees Ubx(i,j) = Vcol_x(i,j)-Ux_deformation(i,j) end if end if coldx end if flotx end do end do do j=1,ny do i=1,nx floty: if (flotmy(i,j)) then coef_defmy(i,j) = 1. Uybar(i,j) = Vcol_y(i,j) ! dmr below is a cast from a real*4 to logical*4 ! dmr cannot be implicit in gfortran ! dmr flgzmy(i,j) = Vcol_y(i,j) flgzmy(i,j) = (nint(Vcol_y(i,j)).ne.0) uydef(i,j) = 0. Uby(i,j) = Vcol_y(i,j) if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143)) call tracebug(' test0') else coldy: if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then !base froide !meme test que dans sliding Bindshadler ddby(i,j) = 0. Uby(i,j) = 0. if (abs(Uy_deformation(i,j)).gt.0.01) then ! calcule par calc_ubar_def coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j) if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143)) call tracebug(' test1') else if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143)) call tracebug(' test2') coef_defmy(i,j)=1. end if else ! base temperee if (abs(Uy_deformation(i,j)).gt.abs(Vcol_y(i,j))) then ! que deformation coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j) ddby(i,j) = 0. Uby(i,j) = 0. if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143)) call tracebug(' test3') else coef_defmy(i,j)=1. ! pas de rescaling de la deformation ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j) ! pb si directions opposees Uby(i,j) = Vcol_y(i,j)-Uy_deformation(i,j) if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.144)) call tracebug(' test4') end if end if coldy end if floty end do end do if (itracebug.eq.1) call tracebug(' apres test floty') ! remarque coef peut être <1 a cause de la pente locale, mais il joue toujours son role debug_3D(:,:,61)=coef_defmx(:,:) debug_3D(:,:,62)=coef_defmy(:,:) calib: do iglen=n1poly,n2poly do k=1,nz SA_mx(:,:,k,iglen) = SA_mx(:,:,k,iglen)*coef_defmx(:,:) S2A_mx(:,:,k,iglen) = S2A_mx(:,:,k,iglen)*coef_defmx(:,:) SA_my(:,:,k,iglen) = SA_my(:,:,k,iglen)*coef_defmy(:,:) S2A_my(:,:,k,iglen) = S2A_my(:,:,k,iglen)*coef_defmy(:,:) end do end do calib end subroutine calc_coef_vitbil !> SUBROUTINE: limit_coef_vitbil !! Calibrate Sa and S2a to force velocity to be balance_velocity in a consistent way. !! @note The difference with calc_coef_vitbil is that the coefficient is limited to the range 0.5-2 !! @note - if > 2, sliding is assumed. !! @note - if < 0.5, it means that ice is actually colder than computed (to be implemented) !! @note Must be called after flowlaw and before velocities !< subroutine limit_coef_vitbil call slope_surf where (abs(sdx(:,:)).lt.1.e-6) sdx(:,:)=-sign(1.e-6,Vcol_x(:,:)) end where where (abs(sdy(:,:)).lt.1.e-6) sdy(:,:)=-sign(1.e-6,Vcol_y(:,:)) end where call calc_ubar_def if (itracebug.eq.1) call tracebug(' Entree dans routine calc_coef_vitbil') !---------------compute coef_defmx ----------------------------------- do j=1,ny do i=1,nx flotx: if (flotmx(i,j)) then coef_defmx(i,j)=1. init_coefmx(i,j)= coef_defmx(i,j) else coldx: if ((coefmxbmelt(i,j).le.0.).and.(.not.gzmx(i,j))) then !base froide !meme test que dans sliding Bindshadler ! sans doute trop fort ddbx(i,j)=0. if (abs(Ux_deformation(i,j)).gt.0.01) then ! calcule par calc_ubar_def coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j) init_coefmx(i,j)= coef_defmx(i,j) if (coef_defmx(i,j).gt.1.) then coef_defmx(i,j)=1. ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j) end if else coef_defmx(i,j)=1. init_coefmx(i,j)= coef_defmx(i,j) end if else ! base temperee if (abs(Ux_deformation(i,j)).gt.abs(Vcol_x(i,j))) then ! only deformation coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j) init_coefmx(i,j)= coef_defmx(i,j) ddbx(i,j) = 0 else coef_defmx(i,j)=1. ! no rescaling of deformation init_coefmx(i,j)= coef_defmx(i,j) ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j) ! pb si directions opposees end if end if coldx end if flotx end do end do do j=1,ny do i=1,nx floty: if (flotmy(i,j)) then coef_defmy(i,j)=1. init_coefmy(i,j)= coef_defmy(i,j) else coldy: if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then !base froide !meme test que dans sliding Bindshadler ! sans doute trop fort ddby(i,j)=0. if (abs(Uy_deformation(i,j)).gt.0.01) then ! calcule par calc_ubar_def coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j) init_coefmy(i,j)= coef_defmy(i,j) if (coef_defmy(i,j).gt.1.) then coef_defmy(i,j)=1. ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j) end if else coef_defmy(i,j)=1. init_coefmy(i,j)= coef_defmy(i,j) end if else ! base temperee if (abs(Uy_deformation(i,j)).gt.abs(Vcol_y(i,j))) then ! only deformation coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j) init_coefmy(i,j)= coef_defmy(i,j) ddby(i,j) = 0 else coef_defmy(i,j)=1. ! no rescaling of deformation init_coefmy(i,j)= coef_defmy(i,j) ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j) ! pb si directions opposees end if end if coldy end if floty end do end do ! remarque coef peut être <1 a cause de la pente locale, mais il joue toujours son role debug_3D(:,:,73)=init_coefmx(:,:) debug_3D(:,:,74)=init_coefmy(:,:) debug_3D(:,:,61)=coef_defmx(:,:) debug_3D(:,:,62)=coef_defmy(:,:) calib: do iglen=n1poly,n2poly do k=1,nz SA_mx(:,:,k,iglen) = SA_mx(:,:,k,iglen)*coef_defmx(:,:) S2A_mx(:,:,k,iglen) = S2A_mx(:,:,k,iglen)*coef_defmx(:,:) SA_my(:,:,k,iglen) = SA_my(:,:,k,iglen)*coef_defmy(:,:) S2A_my(:,:,k,iglen) = S2A_my(:,:,k,iglen)*coef_defmy(:,:) end do end do calib ! call ajust_ghf ATTENTION NE PAS ACTIVER SAUF TESTS ! debug_3D(:,:,73)= init_coefmx(:,:) ! debug_3D(:,:,74)= init_coefmy(:,:) end subroutine limit_coef_vitbil !> SUBROUTINE: calc_ubar_def !! Calculate velocity due to deformation before calibration !< subroutine calc_ubar_def ! calculate velocity due to deformation before calibration implicit none ! extrait de diffusiv pour calculer la partie deformation real :: glenexp real :: inv_4dx ! inverse de dx pour eviter les divisions =1/(4*dx) real :: inv_4dy ! inverse de dy pour eviter les divisions =1/(4*dy) inv_4dx=1./(4*dx) inv_4dy=1./(4*dy) do iglen=n1poly,n2poly glenexp=max(0.,(glen(iglen)-1.)/2.) do j=1,ny do i=1,nx if (.not.flotmy(i,j)) then ddy(i,j,iglen)=((slope2my(i,j)** & ! SLOPE2my calcule dans slope_surf glenexp)*(rog)**glen(iglen)) & *Hmy(i,j)**(glen(iglen)+1) endif if (.not.flotmx(i,j)) then ddx(i,j,iglen)=((slope2mx(i,j)** & ! slope2mx calcule en debut de routine glenexp)*(rog)**glen(iglen)) & *Hmx(i,j)**(glen(iglen)+1) endif end do end do end do do j=2,ny-1 do i=2,nx-1 ux_deformation(i,j)=0. Uy_deformation(i,j)=0. do iglen=n1poly,n2poly ux_deformation(i,j)=ux_deformation(i,j)+ddx(i,j,iglen)*s2a_mx(i,j,1,iglen) Uy_deformation(i,j)=Uy_deformation(i,j)+ddy(i,j,iglen)*s2a_my(i,j,1,iglen) end do end do end do do j=2,ny-1 do i=2,nx-1 ux_deformation(i,j)=ux_deformation(i,j)*(-sdx(i,j)) Uy_deformation(i,j)=Uy_deformation(i,j)*(-sdy(i,j)) end do end do debug_3D(:,:,63)=ux_deformation(:,:) debug_3D(:,:,64)=uy_deformation(:,:) if (itracebug.eq.1) call tracebug(' fin de calc_ubar_def ') end subroutine calc_ubar_def !> SUBROUTINE: ajust_ghf !! Ajuste le flux geothermique pour avoir une temperature coherente !! avec les vitesses de bilan !< subroutine ajust_ghf implicit none real,dimension(nx,ny) :: corr_ghf !< correction du gflux geothermique real,dimension(nx,ny) :: coefdef_maj !< coefficient deformation real :: increment_ghf real :: ghf_min increment_ghf=0.00001 !< exprime en W/m2 ghf_min=0.025 do j=2,ny-1 do i=2,nx-1 if (.not. flot(i,j)) then coefdef_maj(i,j)= (init_coefmx(i,j)+init_coefmx(i+1,j))+ & (init_coefmy(i,j)+init_coefmy(i,j+1)) coefdef_maj(i,j)=0.25*coefdef_maj(i,j) ! un coef trop grand indique eventuellement une base trop froide if ((coefdef_maj(i,j).gt.1.5).and.(ibase(i,j).eq.1)) then ! en base froide ghf(i,j)=ghf(i,j)-SECYEAR*increment_ghf ! attention ghf est negatif ! un coef trop petit indique une base trop chaude else if ((coefdef_maj(i,j).lt.0.7).and.(ghf(i,j).lt.-secyear*ghf_min)) then ghf(i,j)=ghf(i,j)+SECYEAR*increment_ghf ! attention ghf est negatif end if end if end do end do debug_3D(:,:,75)=(ghf0(:,:)-ghf(:,:))*1000./secyear end subroutine ajust_ghf end module spinup_vitbil