!> \file steps_time_loop.f90 !! TOOOOOOOOO DOOOOOOOOOOOOO !! Step_grisli contient !< !-------------------------------------------------------------------------- ! step_time_loop est a l'interieur de la boucle temps. ! il commence par icethick, qui détermine dt, isynchro ! ensuite contient les appels aux diverses sorties ! puis un appel a step_grisli subroutine step_time_loop use module3d_phy use module_choix ! module de choix du type de run ! module_choix donne acces a tous les modules ! de declaration des packages use icetempmod use sorties_ncdf_grisli use flottab_mod use diagno_mod ! use track_debug implicit none integer :: nt_init = 0 ! number of loops for initialisation of ice thickness integer :: iter_visco ! number of iterations for ssa viscosity if (itracebug.eq.1) call tracebug('Entree dans step_time_loop ') if (itracebug.eq.1) write(num_tracebug,*) 'nt = ',nt,' iter_beta = ',iter_beta ! ice thickness evolution !============================= spinup_icethick: if ((ispinup.eq.0.or.ispinup.eq.2.or.nt.lt.nt_init) & .and.(iter_beta.eq.0)) then !------------------------------------------------------------------- ! spinup_icethick : Le calcul changement d'epaisseur (icethick) !---------------- ! - est fait dans tous les cas pour les premieres boucles temporelles ! - pour ispinup = 0 (run standard) sauf si iter_beta = 0 ! - pour ispinup = 2 (force avec vitesses de bilan) ! ! else ! juste un increment pour le temps ! - pour ispinup = 1 et ispinup = 3 temperature calculation ! !------------------------------------------------------------------- if (itracebug.eq.1) call tracebug(' Dans spinup_icethick : appelle icethick') call icethick3() call calving if (itracebug.eq.1) call tracebug('apres calving') ! new calculation of ice surface elevation where(.not.flot(:,:)) S(:,:)=Bsoc(:,:)+H(:,:) B(:,:)=Bsoc(:,:) end where ! calcul de Hmx et Hmy -> shift=-1, dim=1 -> H(i-1,j) hmx(:,:)=0.5*(H(:,:)+eoshift(H(:,:),shift=-1,boundary=0.,dim=1)) hmy(:,:)=0.5*(H(:,:)+eoshift(H(:,:),shift=-1,boundary=0.,dim=2)) hmx(1,:)=0. hmy(:,1)=0. hmx(:,:) = max(hmx(:,:),0.) hmy(:,:) = max(hmy(:,:),0.) ! mise a 0 des vitesses quand epaisseur nulle where (hmx(:,:).le.1.) uxbar(:,:)=0. endwhere where (hmy(:,:).le.1.) uybar(:,:)=0. endwhere else ! time increment when no icethick : ispinup= 1 or 3 ! loop spinup_icethick hdot(:,:)=0. time=nt*dtmax if (itracebug.eq.1) write(num_tracebug,*) ' Dans spinup_icethick: no icethick ',ispinup, & ' time=',time end if spinup_icethick ! end of ice thickness evolution !==================================================================== ! outputs !==================================================================== ! time outputs (volume, extent, ...) !------------------------------------------------------------------ if ((mod(abs(TIME),NDISP*1.).lt.dtmin).or.(isynchro.eq.1).or.(nt.eq.1)) then call shortoutput() endif ! horizontal plan snapshots !------------------------------------------------------------------ ! Premiere partie des sorties horizontales ! if (mod(abs(dble(TIME)),dble(DTSORTIE)).lt.dble(dtmin)) then ! if ((mod(abs(TIME),dtt).lt.dtmin).or.(isynchro.eq.1)) then isynchro=1 ! sorties hz call testsort_time(real(time)) if (iglob_hz.eq.1) then call sortie_hz_multi ! pour les variables déclarées dans 3D ! call hz_output(real(time)) endif else iglob_hz=0 endif ! vertical plan snapshots (profiles) !------------------------------------------------------------------ if ((NT.eq.1) & ! .or.(NT.eq.2) .or.(mod(abs(dble(TIME)),dble(DTPROFILE)).lt.dble(dtmin)) & .or.(abs(TIME-TEND).lt.dtmin)) then ! call sortieprofile() endif ! sorties netcdf Hassine (aout 2010) (2D and 3D) !------------------------------------------------------------------ call testsort_time_ncdf(time) if (iglob_ncdf .EQ. 1) call sortie_ncdf_cat ! sortie compteur tous les dtcpt ans !------------------------------------------------------------------ !iout == 1 sortie cptr !iout == 2 sortie nc if (itracebug.eq.1) call tracebug('dans steps_time_loop avant call out_recovery ') call out_recovery(iout) ! end of outputs !=================================================================== !==================================================================== ! call step_thermomeca() ! ! contient tout le reste de la dynamique et du forçage !==================================================================== end subroutine step_time_loop !--------------------------------------------------------------------- ! The subroutine step_thermomeca() calls the following processes ! ! - climatic forcing and its impact on surface mass balance ! - dynamic velocities (SIA) ! - ice temperature subroutine step_thermomeca() use module3d_phy use module_choix ! module de choix du type de run ! module_choix donne acces a tous les modules ! de declaration des packages use icetempmod use sorties_ncdf_grisli use flottab_mod use diagno_mod use resolmeca_SIA_L1 ! use track_debug implicit none integer :: m logical :: shelf_vitbil = .true. ! si vrai on prend les vitesses de bilan pour le shelf integer :: nt_init_tm = 0 ! number of loops for initialisation of thermo mechanical integer :: iter_visco ! number of iterations for ssa viscosity if (ispinup.le.1) shelf_vitbil = .false. ! general case, ice shelves velocities are computed by diagnoshelf timemax=time if (itracebug.eq.1) write(num_tracebug,*) 'entree dans step_thermomeca', & ' nt=',nt,'iter_beta=',iter_beta !-------------------------------------------------------------------- ! lecture dans un petit fichier "runname.limit" au cas ou les ! parametres du run seraient changes ! !-------------------------------------------------------------------- ! call limit_file(nt,real(time),dt,tend,dtsortie,dtcpt,testdiag,dtt,runname) !if (time.lt.tgrounded) then ! marine=.false. !else marine=.true. !endif ! update the regions (floating, ice ...) !------------------------------------------- call masque() call flottab() if (iter_beta.gt.0) call init_dragging !-----------------------------------------------------------------------! ! debut d'un bloc appels tous les dtt (pas de temps long- asynchrone) ! !------------- -------------------------------------------------------! isync_1: if (ISYNCHRO.eq.1) then ! debut bloc appels tous les dtt ! climatic forcing !===================== call forclim call ablation ! SIA dynamic velocities (deformation + loi de glissement ) !======================== ! remark : velocity is 3D ! to update the velocity for ice temperature calculation ! - surface mass balance for velocity (because ablation was called) ! - impact of surface slope spinup_1_veloc: if ((ispinup.ne.2).or.(nt.lt.nt_init_tm)) then ! cas general !--------------------------------------------------------------------------------- ! spinup_1 : calcul des vitesses dynamiques : !------------ ! - est fait dans tous les cas pour les premieres boucles temporelles ! - est fait pour ispinup=3 pour calibration forme de vitesse ! - n'est pas fait pour ispinup=2 (test balance vel sur icethick) ! !--------------------------------------------------------------------------------- if (itracebug.eq.1) call tracebug(' Dans spinup_1_veloc : calculer SIA_velocities') call Neffect() call SIA_velocities() else ! ispinup = 2 call force_balance_vel ! fait l'ensemble du champ (flgz ou pas) endif spinup_1_veloc ! ice temperature !=================================================================== ! update values in the structures Geom_g, Temperature_g, ... calc_temp: if ((nt.gt.2).and.(geoplace(1:5).ne.'mism3')) then if (itracebug.eq.1) call tracebug('avant appel icetemp') call icetemp ! subroutines pour le calcul de la fusion basale call bmeltshelf call bmelt_grounded endif calc_temp ! flowlaw !=================================================================== spinup_2_flowlaw: if ((ispinup.ne.2).or.(nt.lt.nt_init_tm)) then ! cas general !--------------------------------------------------------------------------------- ! spinup_2_flowlaw !---------------- ! si ispinup = 2 , on ne fait pas le calcul de flowlaw, ni de diffusiv ! puisque les vitesses de bilan sont imposees ! (initialisee nt < nt_init_tm) ! ! Flowlaw et diffusiv : servent soit directement (ispinup = 0,1) soit pour le ! rescaling (ispinup = 3) !--------------------------------------------------------------------------------- if (itracebug.eq.1) call tracebug(' Dans spinup_2_flowlaw') call flow_general do iglen=n1poly,n2poly call flowlaw(iglen) end do end if spinup_2_flowlaw if (ispinup.eq.3) then call calc_coef_vitbil ! recalibre sa_mx et s2a_mx et ddbx pour toute la zone posee end if ! Dans le shelf flgz = Ubar = Vbil et udef=0 ! isostasy !================= spinup_3_bed: if (ispinup.eq.0.or.nt.lt.nt_init_tm) then !--------------------------------------------------------------------------------- ! spinup_3_bed !---------------- ! Ne passe dans ce bloc que quand on veut l'isostasie, donc si variations epaisseur ! run standard (ispinup=0) et pour les initialisations !--------------------------------------------------------------------------------- if (itracebug.eq.1) call tracebug(' Dans spinup_3_bed') if ((nbed.eq.1).and.nt.ne.1.and.isynchro.eq.1) then call bedrock() ! bedrock adjustment endif call lake_to_slv ! pour traiter dynamiquement les shelves sur lac end if spinup_3_bed call tracer ! aurel neem (attention, position a verifier) endif isync_1 ! fin du bloc tous les dtt isync_1 !------------------------------------ spinup_4_vitdyn: if ((ispinup.le.1).or.(nt.lt.nt_init_tm)) then !--------------------------------------------------------------------------------- ! spinup_4_vitdyn !---------------- ! Passe dans ce bloc que quand on veut calculer les vitesses dynamiques : ! - initialisations (nt.lt.nt_init_tm) ! - run standard (ispinup=0) (y compris iter_beta.ne.0) ! - spinup temperature avec vitesses dynamiques (ispinup = 1) !--------------------------------------------------------------------------------- call Neffect() call diffusiv() iterbeta: if (iter_beta.ne.0) then ! only for iterations on beta if (.not.shelf_vitbil) then ! bloc si les vitesses shelves viennent de diagnoshelf ! Vcol declare dans le module spinup_mod ! ou dans le dragging si no_spinup where (flotmx(:,:)) uxbar(:,:)=uxflgz(:,:) elsewhere uxbar(:,:)=Vcol_x(:,:) end where where (flotmy(:,:)) uybar(:,:)=uyflgz(:,:) elsewhere uybar(:,:)=Vcol_y(:,:) end where else uxbar(:,:)=Vcol_x(:,:) uybar(:,:)=Vcol_y(:,:) end if end if iterbeta ! strain rate calculations !=========================== call strain_rate() if (iter_beta.ne.0) then ! on n'a utilise les vitesses de bilan que pour strain_rate "Dirichlet ?" uxbar(:,:)=uxflgz(:,:)+uxdef(:,:) uybar(:,:)=uyflgz(:,:)+uydef(:,:) end if ! ice shelves and ice streams velocities : SSA !============================================== !----------------------------------------------------------------------- ! debut d'un deuxieme bloc appels tous les dtt !------------- ------------------------------------------------------- isync_2: if ((shelfy).and.(marine).and.(isynchro.eq.1)) then !les 3 lignes ci-dessous pour avoir les champs avant diagno (en cas de pb colineaire) !call init_sortie_ncdf !call sortie_ncdf_cat if (icompteur.eq.1) then ! reprise de tout le vecteur d'etat iter_visco = 1 else if ((icompteur.ne.1).and.(nt.le.1)) then ! initialisation avec reprise partielle ! ou nulle du vecteur d' etat iter_visco = 10 else iter_visco = 2 end if do m=1,iter_visco call diagnoshelf call mix_SIA_L1 call strain_rate() end do ! iterations pour affiner beta en fonction des vitesses de bilan if ((time.ge.time_iter).and.(time.le.time_iter+15.)) then !iterations seulment pendant 15 ans do m=1,nb_iter_vitbil call beta_iter_vitbil call diagnoshelf call mix_SIA_L1 call strain_rate() end do end if ! fin des iterations beta_iter_vitbil endif isync_2 call mix_SIA_L1 end if spinup_4_vitdyn end subroutine step_thermomeca