!> \file calving_frange.f90 !! Module et routines qui calculent le calving avec la methode de Vincent !! Calving adapte pour simulation abuk ABUMIP : suppression de toute la glace !! flottante via le calving a chaque pas de temps. !< !> \namespace module3D_phy !! Module et routines qui calculent le calving avec la methode de Vincent !! \author ... !! \date ... !! @note recopie direct de icethick5-ant. A revoir en particulier, accroche avec la glace posee !! @note Used modules !! @note - use module3D_phy !! @todo il faudrait que le calving soit inactif dans les zones ou l epaisseur est imposee !< module calving_frange_abuk use module3D_phy use bilan_eau_mod implicit none real, dimension (nx,ny) :: hmhc ! hauteur au dessus de la coupure real, dimension (nx,ny) :: bil_tot ! bilan surface et fond (bm-bmelt) real, dimension (nx,ny) :: hcoup ! epaiseur de coupure au temps time real :: hcoup_plateau ! coupure points peu profonds real :: hcoup_abysses ! coupure points ocean profond real :: prof_plateau ! profondeur max des points peu profonds real :: prof_abysses ! profondeur min des points ocean profond integer :: meth_hcoup ! pour avoir hcoup dépendant du coefbmshelf integer :: ifrange ! 0 pas de traitement particulier pres du bord, 1 -> franges integer :: iin2,jin2,iin3,jin3 ! pour la detection polynies logical :: testmij,testpij,testimj,testipj logical :: bilan_surf_fond ! vrai si bm-bmelt est positif logical :: avalw,avale,avals,avaln,interieur contains !--------------------------------------------------------------------------------------- subroutine init_calving namelist/calving/Hcoup_plateau,Hcoup_abysses,prof_plateau,prof_abysses,ifrange,meth_hcoup calv(:,:)=0. calv_dtt(:,:)=0. ! formats pour les ecritures dans 42 428 format(A) ! lecture des parametres du run block eaubasale1 !-------------------------------------------------------------------- rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,calving) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '&calving ! nom du bloc calving méthode Vincent ' write(num_rep_42,*) write(num_rep_42,*) 'Hcoup_plateau = ', Hcoup_plateau write(num_rep_42,*) 'Hcoup_abysses = ', Hcoup_abysses write(num_rep_42,*) 'prof_plateau = ', prof_plateau write(num_rep_42,*) 'prof_abysses = ', prof_abysses write(num_rep_42,*) 'ifrange = ', ifrange write(num_rep_42,*) 'meth_hcoup = ', meth_hcoup write(num_rep_42,*)'/' write(num_rep_42,428) '! Hcoup epaisseurs de coupure pour les zones peu prodondes et profondes' write(num_rep_42,428) '! Hcoup_plateau pas de traitement particulier sur les bords' write(num_rep_42,428) '! ifrange=1 -> traitement de Vincent avec ice shelves frangeants partout' write(num_rep_42,428) '! ifrange=2 -> ice shelves frangeants seulement si bm-bmelt positif' write(num_rep_42,*) '! meth_hcoup pour faire eventuellement varier Hcoup avec le climat' write(num_rep_42,*) '! Hcoup_plateau=',Hcoup_plateau write(num_rep_42,*) '! Hcoup_abysses=',Hcoup_abysses write(num_rep_42,*) '! prof_plateau=',prof_plateau write(num_rep_42,*) '! prof_abysses=',prof_abysses write(num_rep_42,*) ! afq -- coupure depend de la profondeur: ! ! hcoup prof_abysses ! ^ v ! | _______ hcoup_abysses ! | / ! | / ! | / ! | / ! | hcoup_plateau _______/ ! ^ ! prof_plateau Hcoup(:,:) = min ( max( & (-(Bsoc0(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau) & *(hcoup_abysses-hcoup_plateau)+hcoup_plateau & , hcoup_plateau), hcoup_abysses ) if (meth_hcoup.eq.1) then Hcoup(:,:)=coefbmshelf*Hcoup(:,:) Hcoup(:,:)=min( max(Hcoup (:,:),Hcoup_plateau),Hcoup_abysses) else if (meth_hcoup.eq.2) then Hcoup(:,:)=coefbmshelf*Hcoup(:,:) Hcoup(:,:)=max(Hcoup(:,:),Hcoup_plateau) else if (meth_hcoup.eq.3) then Hcoup(:,:)=coefbmshelf*Hcoup(:,:) Hcoup(:,:)=min(max(Hcoup(:,:),0.),Hcoup_abysses) endif end subroutine init_calving !--------------------------------------------------------------------------------------- subroutine calving integer :: I_did_something ! pour la boucle sur le calving ! initialisation calving calv(:,:)=0. ! calcul du dhdt lagrangien dhdt(:,:)=bm(:,:)-bmelt(:,:)-H(:,:)*(epsxx(:,:)+epsyy(:,:)) ! calcul du bilan surface et fond : divise par 2 car utilise dans des moyennes bil_tot(:,:)=0.5*(bm(:,:)-bmelt(:,:)) Hcoup(:,:) = min ( max( & (-(Bsoc(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau) & *(hcoup_abysses-hcoup_plateau)+hcoup_plateau & , hcoup_plateau), hcoup_abysses ) if (meth_hcoup.eq.1) then Hcoup(:,:)=coefbmshelf*Hcoup(:,:) Hcoup(:,:)=min( max(Hcoup (:,:),Hcoup_plateau),Hcoup_abysses) else if (meth_hcoup.eq.2) then Hcoup(:,:)=coefbmshelf*Hcoup(:,:) Hcoup(:,:)=max(Hcoup(:,:),Hcoup_plateau) else if (meth_hcoup.eq.3) then Hcoup(:,:)=coefbmshelf*Hcoup(:,:) Hcoup(:,:)=min(max(Hcoup(:,:),0.),Hcoup_abysses) endif ! hauteur au dessus de la coupure hmhc(:,:)=H(:,:)-Hcoup(:,:) ! coupure de l'ice shelf !--------------------------- ! on divise le test en 2 'ifext' et 'ifint' pour reduire les tests redondants MULTI_CALV_LOOP : do l=1,10 ! calcul du dhdt lagrangien dhdt(:,:)=bm(:,:)-bmelt(:,:)-H(:,:)*(epsxx(:,:)+epsyy(:,:)) ! hauteur au dessus de la coupure hmhc(:,:)=H(:,:)-Hcoup(:,:) I_did_something = 0 do j=2,ny-1 do i=2,nx-1 ifext: if((flot(i,j)).and.(h(i,j).le.hcoup(i,j)).and.(h(i,j).gt.0.)) then ! ifext: pour les noeuds flottants englaces avec h < hcoup !ifint: if((front(i,j).gt.0).and.(front(i,j).lt.4)) then !cdc pb avec front ifint: if((front(i,j).lt.4).or.((front(i-1,j)+front(i+1,j)+front(i,j-1)+front(i,j+1)).lt.16)) then ifint: if((H(i-1,j).lt.2.).or.(H(i+1,j).lt.2.).or.(H(i,j-1).lt.2.).or.(H(i,j+1).lt.2)) then ! si on est au bord avec test sur H ! ifint: le point doit avoir au - 1 voisin mais ne pas etre entouré de glace ! ce qui evite la formation des polynies dans les shelfs ! on regarde dhdt minimum pour voir si 1 des voisins peut alimenter le point ! hmhc est l'épaisseur en plus de l'épaisseur de coupure ! ux/dx=1/(deltat) ou deltat est le temps nécessaire ! à la glace pour passer d'un noeud à l'autre ! la deuxieme partie du test revient à dh/dt*deltat > hmhc ! Rajout vince : si le point amont est pose -> test vrai. ! rappel des differents ifrange !------------------------------- !Attention ifrange 1,2 ont sans doute un bug (il faudrait abs(uxbar)) ! ifrange=1 Rajout vince : si un point voisin est pose -> test vrai. ! comme dans article fenno2007 ! ifrange=2 pour les points ayant des voisins poses, garde le point ! seulement si le bilan surface fond est positif. ! bm(i,j)-bmelt(i,j) > 0. ! ifrange=3 pas de test sur l'epaisseur du point amont. ! test dhdt > -hmhc/deltat= -hmhc abs(U)/dx (correction du bug) ! + test sur bm(i,j)-bmelt(i,j) > 0. ! ifrange = 4 idem 3 mais avec test sur l'epaisseur du point amont ! ifrange = 0 pas de traitement specifique près du continent ! ifrange = -1 ancienne version MIS11 ! ifrange = 5 semble idem 0 ! ifrange = 6 ?? ! ifrange = 7 calving drastique sauf si coefbmshelf < 0.5 if (ifrange.eq.1) then ! Rajout vince : si un point voisin est pose -> test vrai. ! comme dans article fenno2007 testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.) & ! voisin (i-1,j) amont et > hcoup .and.(dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx)))& .or.(.not.flot(i-1,j)) ) ! testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx))) & .or.(.not.flot(i+1,j)) ) ! testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.) & ! voisin (i,j-1) amont et > hcoup .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx)))& .or.(.not.flot(i,j-1)) ) ! testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx)))& .or.(.not.flot(i,j+1)) ) ! else if (ifrange.eq.2) then ! pour les points ayant des voisins posés, seulement si le bilan ! surface fond est positif. bm(i,j)-bmelt(i,j) > 0. bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.) testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.) & ! voisin (i-1,j) amont et > hcoup .and. (dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx))) & .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) ! testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx))) & .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) ! testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.) & ! voisin (i,j-1) amont et > hcoup .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx)))& .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) ! testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx)))& .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) ! else if (ifrange.eq.3) then ! nouvelle formulation Cat mars 08 ! pas de test sur l'epaisseur du point amont. ! test dhdt > -hmhc/deltat= -hmhc abs(U)/dx ! pour les points ayant des voisins posés, seulement si le bilan ! surface fond est positif. bm(i,j)-bmelt(i,j) > 0. bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.) testmij=( (uxbar(i,j).ge.0.) & ! voisin (i-1,j) amont et > hcoup .and.(dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j))/dx))) & .or.((.not.flot(i-1,j)).and.bilan_surf_fond ) ! testpij=( (uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx))) & .or.((.not.flot(i+1,j)).and.bilan_surf_fond ) ! testimj=(( uybar(i,j).ge.0.) & ! voisin (i,j-1) amont et > hcoup .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx))) & .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ! testipj=((uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx))) & .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ! else if (ifrange.eq.4) then ! idem 3 mais avec test sur l'épaisseur du point amont bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.) testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.) & ! voisin (i-1,j) amont et > hcoup .and. (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) & .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) ! testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) & .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) ! testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.) & ! voisin (i,j-1) amont et > hcoup .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))& .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) ! testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))& .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) ! else if (ifrange.eq.0) then ! pas de traitement special pres du continent testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.) & ! voisin (i-1,j) amont et > hcoup .and.(dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx))) testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.) & ! voisin (i,j-1) amont et > hcoup .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx))) testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx))) else if (ifrange.eq.-1) then ! ancienne version MIS11 testmij=((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.) & ! voisin (i-1,j) amont et > hcoup .and.(dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx))) testpij=((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx))) testimj=((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.) & ! voisin (i,j-1) amont et > hcoup .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx))) testipj=((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx))) else if (ifrange.eq.5) then ! pas de traitement special pres du continent testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.) & ! voisin (i-1,j) amont et > hcoup .and.((bil_tot(i-1,j)+bil_tot(i,j)) & .gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup .and.((bil_tot(i+1,j)+bil_tot(i,j)) & .gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx))) testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.) & ! voisin (i,j-1) amont et > hcoup .and.((bil_tot(i,j-1)+bil_tot(i,j)) & .gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx))) testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup .and.((bil_tot(i,j+1)+bil_tot(i,j)) & .gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx))) else if (ifrange.eq.6) then ! pas de traitement special pres du continent testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.) & ! voisin (i-1,j) amont et > hcoup .and.(2.*bil_tot(i,j) & .gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup .and.(2.*bil_tot(i,j) & .gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx))) testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.) & ! voisin (i,j-1) amont et > hcoup .and.(2.*bil_tot(i,j) & .gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx))) testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup .and.(2.*bil_tot(i,j) & .gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx))) else if (ifrange.eq.7) then ! drastique calving sauf si coefbmelt.lt.0.5 testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.) & ! voisin (i-1,j) amont et > hcoup .and.(coefbmshelf.lt.1.)) testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup .and.(coefbmshelf.lt.1.)) testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.) & ! voisin (i,j-1) amont et > hcoup .and.(coefbmshelf.lt.1.)) testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup .and.(coefbmshelf.lt.1.)) endif ! detection des polynies dans les shelfs. on regarde vers l'aval ! 1 des 3 voisins aval > hcoup du cote west (i-1) iin2=max(1,i-2) iin3=max(1,i-3) avalw=((hmhc(i-1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) & .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i,j).le.0.) ! 1 des 3 voisins aval > hcoup du cote est (i+1) iin2=min(i+2,nx) iin3=min(i+3,nx) avale=((hmhc(i+1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) & .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i+1,j).ge.0) ! 1 des 3 voisins aval > hcoup du cote sud (j-1) jin2=max(1,j-2) jin3=max(1,j-3) avals=((hmhc(i,j-1).gt.0.).or.(hmhc(i,jin2).gt.0.) & .or.(hmhc(i,jin3).gt.0.)).and.(uybar(i,j).le.0.) ! 1 des 3 voisins aval > hcoup du cote nord (j-1) jin2=min(j+2,ny) jin3=min(j+3,ny) avaln=((hmhc(i,j+1).gt.0.).or.(hmhc(i,jin2).gt.0.) & .or.(hmhc(i,jin2).gt.0.)) .and.(uybar(i,j+1).ge.0.) interieur=(avalw.or.avale).and.(avals.or.avaln) if ((.not.(testmij.or.testpij.or.testimj.or.testipj)) & ! pas suffisament alimente .and.(.not.interieur)) then ! et pas interieur calv(i,j)=-h(i,j) !cdc H(i,j)=1. !cdc 1m H(i,j)=min(1.,max(0.,(sealevel - Bsoc(i,j))*row/ro-0.01)) H(i,j)=0. S(i,j)=H(i,j)*(1.-ro/row) + sealevel B(i,j)=S(i,j) - H(i,j) ! ATTENTION ne pas mettre ice=0 sinon degradation bilan d'eau (bm et bmelt non comptabilises dans ce cas) I_did_something = I_did_something + 1 ! if (l.ge.2) then ! print*,'calving l ij',l,i,j ! endif endif end if ifint end if ifext end do end do if (I_did_something.eq.0) then ! print*,'stop MULTI_CALV_LOOP l=',l EXIT MULTI_CALV_LOOP else ! print*,'calving continue l',l,I_did_something endif enddo MULTI_CALV_LOOP ! on met en calving les points detectes iceberg : !~ where (iceberg(:,:).and.(H(:,:).gt.0.)) !~ calv(:,:)=-h(:,:) !~ ice(:,:)=0 !~ H(:,:)=0. !~ S(:,:)=H(:,:)*(1.-ro/row) + sealevel !~ B(:,:)=S(:,:) - H(:,:) !~ endwhere ! version abumip abuk (kill ice-shelves) where (flot(:,:).and.(H(:,:).gt.0.)) calv(:,:)=-h(:,:) ice(:,:)=0 H(:,:)=0. S(:,:)=H(:,:)*(1.-ro/row) + sealevel B(:,:)=S(:,:) - H(:,:) endwhere calv_dtt(:,:) = calv_dtt(:,:) + calv(:,:) ! somme du calving sur dtt ! calv_dtt est remis a 0 dans steps_time_loop (tous les dtt) end subroutine calving !------------------------------------------------------------------------------------------ end module calving_frange_abuk