source: branches/iLoveclim/SOURCES/steps_time_loop.f90 @ 123

Last change on this file since 123 was 123, checked in by aquiquet, 7 years ago

Merged branch iLOVECLIM to trunk at rev 121

File size: 14.5 KB
Line 
1!> \file steps_time_loop.f90
2!! TOOOOOOOOO DOOOOOOOOOOOOO
3!! Step_grisli contient
4!<
5
6!--------------------------------------------------------------------------
7! step_time_loop est a l'interieur de la boucle temps.
8! il commence par icethick, qui détermine dt, isynchro
9
10! ensuite contient les appels aux diverses sorties
11! puis un appel a step_grisli
12
13
14subroutine step_time_loop       
15
16  use module3d_phy
17  use module_choix ! module de choix du type de run
18                   !  module_choix donne acces a tous les modules
19                   !  de declaration des packages
20  use sorties_ncdf_grisli
21  use flottab_mod
22  use icetempmod
23  use diagno_mod
24  use bilan_eau_mod 
25  !  use track_debug
26
27  implicit none
28  integer :: nt_init = 0      ! number of loops for initialisation of ice thickness
29  integer :: iter_visco       ! number of iterations for ssa viscosity
30
31  if (itracebug.eq.1)  call tracebug('Entree dans step_time_loop ')
32  if (itracebug.eq.1) write(num_tracebug,*) 'nt = ',nt,'  iter_beta = ',iter_beta
33
34
35  ! ice thickness evolution
36  !=============================
37
38
39  spinup_icethick: if ((ispinup.eq.0.or.ispinup.eq.2.or.nt.lt.nt_init) &
40       .and.(iter_beta.eq.0)) then 
41
42     !-------------------------------------------------------------------
43     ! spinup_icethick : Le calcul changement d'epaisseur (icethick)
44     !----------------
45     !   - est fait dans tous les cas pour les premieres boucles temporelles
46     !   - pour ispinup = 0 (run standard) sauf si iter_beta = 0
47     !   - pour ispinup = 2 (force avec vitesses de bilan)
48     !
49     ! else
50     !   juste un increment pour le temps
51     !    -  pour ispinup = 1 et  ispinup = 3  temperature calculation
52     
53     !-------------------------------------------------------------------
54
55     if (itracebug.eq.1)  call tracebug(' Dans spinup_icethick : appelle icethick')
56
57     call icethick3
58
59     call flottab
60
61     call calving
62
63     call ablation_bord
64
65     call bilan_eau
66     if (isynchro.eq.1) then
67         call shortoutput
68         diff_H = 0.
69         Bm_dtt(:,:) = 0.
70         bmelt_dtt(:,:) = 0.
71         calv_dtt(:,:)=0.
72         ablbord_dtt(:,:)=0.
73         diff_H_2D(:,:)=0.
74      endif
75     
76     call flottab
77
78     if (itracebug.eq.1)  call tracebug('apres calving')
79
80     ! new calculation of ice surface elevation
81     where(.not.flot(:,:))
82        S(:,:)=Bsoc(:,:)+H(:,:)
83        S(:,:)=max(S(:,:),sealevel)
84        B(:,:)=Bsoc(:,:)
85     end where
86
87     ! calcul de Hmx et Hmy -> shift=-1, dim=1 -> H(i-1,j)
88
89     hmx(:,:)=0.5*(H(:,:)+eoshift(H(:,:),shift=-1,boundary=0.,dim=1))
90     hmy(:,:)=0.5*(H(:,:)+eoshift(H(:,:),shift=-1,boundary=0.,dim=2))
91     hmx(1,:)=0.
92     hmy(:,1)=0.
93
94     ! mise a 0 des vitesses quand epaisseur nulle
95     where (hmx(:,:).le.1.)
96        uxbar(:,:)=0.
97     endwhere
98     where (hmy(:,:).le.1.)
99        uybar(:,:)=0.
100     endwhere
101
102  else   ! time increment when no icethick : ispinup= 1 or 3    ! loop spinup_icethick
103
104     hdot(:,:)=0.
105     time=nt*dtmax
106
107     if (itracebug.eq.1) write(num_tracebug,*) ' Dans spinup_icethick:  no icethick ',ispinup, &
108          '    time=',time
109
110
111  end if spinup_icethick
112  ! end of ice thickness evolution
113  !====================================================================
114
115
116
117  ! outputs
118  !====================================================================
119
120  ! horizontal plan snapshots
121  !------------------------------------------------------------------
122  !  Premiere partie des sorties horizontales
123
124  !     if (mod(abs(dble(TIME)),dble(DTSORTIE)).lt.dble(dtmin)) then
125  !
126
127
128  if ((mod(abs(TIME),dtt).lt.dtmin).or.(isynchro.eq.1)) then
129     isynchro=1
130
131     !  sorties hz
132     call testsort_time(real(time))
133     if (iglob_hz.eq.1) then
134        call sortie_hz_multi   ! pour les variables déclarées dans 3D
135        !              call hz_output(real(time))
136     endif
137
138  else
139     iglob_hz=0
140
141  endif
142
143
144  !   vertical plan snapshots (profiles)
145  !------------------------------------------------------------------
146
147  !if ((NT.eq.1) &    !   .or.(NT.eq.2)
148  !     .or.(mod(abs(dble(TIME)),dble(DTPROFILE)).lt.dble(dtmin)) &
149  !     .or.(abs(TIME-TEND).lt.dtmin)) then
150
151     !        call sortieprofile()
152
153  !endif
154
155  !  sorties netcdf Hassine (aout 2010)  (2D and 3D)
156  !------------------------------------------------------------------
157
158  call testsort_time_ncdf(time)
159  if (iglob_ncdf .EQ. 1) call sortie_ncdf_cat
160
161
162  !   sortie compteur tous les dtcpt ans
163  !------------------------------------------------------------------
164  !iout == 1 sortie cptr
165  !iout == 2 sortie nc
166  if (itracebug.eq.1)  call tracebug('dans steps_time_loop avant call out_recovery ')
167
168  call out_recovery(iout)
169 
170  ! end of outputs
171
172  !====================================================================
173  !
174  call step_thermomeca()
175  !
176  ! contient tout le reste de la dynamique et du forçage
177  !====================================================================
178
179end subroutine step_time_loop
180
181
182
183
184
185
186!---------------------------------------------------------------------
187! The subroutine  step_thermomeca() calls the following processes
188!
189! - climatic forcing and its impact on surface mass balance
190! - dynamic velocities (SIA)
191! - ice temperature
192
193subroutine step_thermomeca()
194
195
196  use module3d_phy
197  use module_choix ! module de choix du type de run
198  !  module_choix donne acces a tous les modules
199  !  de declaration des packages
200  use icetempmod
201  use sorties_ncdf_grisli
202  use flottab_mod
203  use diagno_mod
204  use resolmeca_SIA_L1
205!  use track_debug
206
207
208  implicit none
209
210  integer :: m
211  logical :: shelf_vitbil = .true.   ! si vrai on prend les vitesses de bilan pour le shelf
212  integer :: nt_init_tm = 0          ! number of loops for initialisation of thermo mechanical
213  integer :: iter_visco              ! number of iterations for ssa viscosity
214
215  if (ispinup.le.1) shelf_vitbil = .false.         ! general case, ice shelves velocities are computed by diagnoshelf
216
217  timemax=time
218  if (itracebug.eq.1) write(num_tracebug,*) 'entree dans step_thermomeca', &
219       '   nt=',nt,'iter_beta=',iter_beta
220
221  !--------------------------------------------------------------------
222  !     lecture dans un petit fichier "runname.limit" au cas ou les
223  !     parametres du run seraient changes 
224  !                                                   
225  !--------------------------------------------------------------------
226
227!  call limit_file(nt,real(time),dt,tend,dtsortie,dtcpt,testdiag,dtt,runname)
228
229  !if (time.lt.tgrounded) then
230  !   marine=.false.
231  !else
232  marine=.true.
233  !endif
234
235  ! update the regions (floating, ice ...)
236  !-------------------------------------------
237  call masque()
238
239  if (iter_beta.gt.0)   call init_dragging
240
241
242  !-----------------------------------------------------------------------!
243  !   debut d'un bloc appels tous les dtt (pas de temps long- asynchrone) !
244  !-------------   -------------------------------------------------------!
245
246
247
248  isync_1: if (ISYNCHRO.eq.1) then          ! debut bloc appels tous les dtt   
249
250     ! climatic forcing
251     !=====================
252     call forclim                             
253     call ablation
254
255     ! oceanic forcing
256     !=====================
257     call bmeltshelf  ! afq -- icetemp is not called before the 4th dtt, so a first call to bmelt here
258
259     ! SIA dynamic velocities  (deformation + loi de glissement )
260     !========================
261
262     ! remark : velocity is 3D
263     !          to update the velocity for ice temperature calculation
264     !            - surface mass balance for velocity (because ablation was called)
265     !            - impact of surface slope
266
267
268     spinup_1_veloc: if ((ispinup.ne.2).or.(nt.lt.nt_init_tm)) then        ! cas general
269
270        !---------------------------------------------------------------------------------
271        ! spinup_1 : calcul des vitesses dynamiques :
272        !------------
273        !            - est fait dans tous les cas pour les premieres boucles temporelles
274        !            - est fait pour ispinup=3 pour calibration forme de vitesse
275        !            - n'est pas fait pour ispinup=2 (test balance vel sur icethick)
276        !
277        !---------------------------------------------------------------------------------
278
279        if (itracebug.eq.1)  call tracebug(' Dans spinup_1_veloc : calculer SIA_velocities')
280
281        call Neffect()                                             
282        call SIA_velocities()
283
284     else    ! ispinup = 2
285
286        call force_balance_vel   ! fait l'ensemble du champ (flgz ou pas)
287
288     endif spinup_1_veloc
289
290
291     ! ice temperature
292     !===================================================================
293
294     ! update values in the structures Geom_g, Temperature_g, ...
295
296     calc_temp:     if ((nt.gt.2).and.(geoplace(1:5).ne.'mism3')) then
297        if (itracebug.eq.1)  call tracebug('avant appel icetemp')
298        call icetemp
299        ! subroutines pour le calcul de la fusion basale
300        ! afq -- 28-10-2016: there is a reset of bmelt in temp_col in the
301        ! current version of the code, so shelf and grounded melt have to
302        ! be called here. TODO: remove reset in icetemp routines??
303        call bmeltshelf
304        call bmelt_grounded
305
306     endif calc_temp
307
308
309     ! flowlaw
310     !===================================================================
311
312     spinup_2_flowlaw: if ((ispinup.ne.2).or.(nt.lt.nt_init_tm)) then    ! cas general
313
314        !---------------------------------------------------------------------------------
315        ! spinup_2_flowlaw
316        !----------------
317        ! si ispinup =   2 , on ne fait pas le calcul de flowlaw, ni de diffusiv
318        !                    puisque les vitesses de bilan sont imposees
319        !                                             (initialisee nt < nt_init_tm)
320        !                   
321        ! Flowlaw et diffusiv : servent soit directement (ispinup = 0,1) soit pour le
322        !                       rescaling (ispinup = 3)
323        !---------------------------------------------------------------------------------
324
325        if (itracebug.eq.1)  call tracebug(' Dans spinup_2_flowlaw')
326
327           call flow_general
328           do iglen=n1poly,n2poly
329              call flowlaw(iglen)
330           end do
331
332
333     end if spinup_2_flowlaw
334
335     if (ispinup.eq.3) then
336        call calc_coef_vitbil                    ! recalibre sa_mx et s2a_mx et ddbx pour toute la zone posee
337     end if                                      ! Dans le shelf flgz = Ubar = Vbil et udef=0
338
339
340     ! isostasy
341     !=================
342
343
344     spinup_3_bed: if (ispinup.eq.0.or.nt.lt.nt_init_tm) then
345        !---------------------------------------------------------------------------------
346        ! spinup_3_bed
347        !----------------
348        ! Ne passe dans ce bloc que quand on veut l'isostasie, donc si variations epaisseur
349        ! run standard  (ispinup=0) et pour les initialisations                 
350        !---------------------------------------------------------------------------------
351
352        if (itracebug.eq.1)  call tracebug(' Dans spinup_3_bed')
353
354!        if ((nbed.eq.1).and.nt.ne.1.and.isynchro.eq.1) then
355         if ((nbed.eq.1).and.nt.ne.1.and.isynchro.eq.1.and.(mod(abs(TIME),50.).lt.dtmin)) then
356           call bedrock                                              !  bedrock adjustment
357        endif
358
359        call lake_to_slv                                             ! pour traiter dynamiquement les shelves sur lac
360
361     end if spinup_3_bed
362
363     call tracer      ! aurel neem (attention, position a verifier)   
364
365  endif isync_1
366
367
368  ! fin du bloc tous les dtt isync_1
369  !------------------------------------
370
371
372  spinup_4_vitdyn: if  ((ispinup.le.1).or.(nt.lt.nt_init_tm)) then 
373
374     !---------------------------------------------------------------------------------
375     ! spinup_4_vitdyn
376     !----------------
377     ! Passe dans ce bloc que quand on veut calculer les vitesses dynamiques :
378     ! - initialisations (nt.lt.nt_init_tm)
379     ! - run standard  (ispinup=0) (y compris iter_beta.ne.0)
380     ! - spinup temperature avec vitesses dynamiques  (ispinup = 1)               
381     !---------------------------------------------------------------------------------
382
383
384
385     call Neffect()
386
387     call diffusiv()
388
389
390     iterbeta:     if (iter_beta.ne.0) then          ! only for iterations on beta
391
392        if (.not.shelf_vitbil) then    !  bloc si les vitesses shelves viennent de diagnoshelf
393
394           !  Vcol declare dans le module spinup_mod
395           !  ou dans le dragging si no_spinup
396
397           where (flotmx(:,:))
398              uxbar(:,:)=uxflgz(:,:)
399           elsewhere
400              uxbar(:,:)=Vcol_x(:,:)
401           end where
402
403           where (flotmy(:,:))
404              uybar(:,:)=uyflgz(:,:)
405           elsewhere
406              uybar(:,:)=Vcol_y(:,:)
407           end where
408        else
409           uxbar(:,:)=Vcol_x(:,:)
410           uybar(:,:)=Vcol_y(:,:)
411        end if
412     end if iterbeta
413
414
415     ! strain rate calculations
416     !===========================
417
418     call strain_rate()
419
420
421     if (iter_beta.ne.0) then               ! on n'a utilise les vitesses de bilan que pour strain_rate "Dirichlet ?"
422        uxbar(:,:)=uxflgz(:,:)+uxdef(:,:)
423        uybar(:,:)=uyflgz(:,:)+uydef(:,:)
424     end if
425
426
427
428     ! ice shelves and ice streams velocities : SSA
429     !==============================================
430
431     !-----------------------------------------------------------------------
432     !              debut d'un deuxieme bloc appels tous les dtt                         
433     !-------------   -------------------------------------------------------
434
435     isync_2:    if ((shelfy).and.(marine).and.(isynchro.eq.1)) then               
436
437        !les 3 lignes ci-dessous pour avoir les champs avant diagno (en cas de pb colineaire)
438        !call init_sortie_ncdf
439        !call sortie_ncdf_cat
440
441        if (icompteur.eq.1) then    ! reprise de tout le vecteur d'etat   
442           iter_visco = 1
443
444        else if ((icompteur.ne.1).and.(nt.le.1)) then  ! initialisation avec reprise partielle
445                                                       ! ou nulle du vecteur d' etat
446          iter_visco =  10
447        else
448           iter_visco = 2
449        end if
450 
451        iter_visco= 10                                 ! warning, test sur l'impact iteration dragging et visco
452                                                       ! Cat 18 janv 2013
453
454           do m=1,iter_visco             
455              call diagnoshelf
456              call  mix_SIA_L1 
457              call strain_rate()
458              if (iter_beta.eq.0) then
459                call dragging         
460              endif
461              if (test_iter_diag.lt.1.e-2) exit
462           end do
463        endif isync_2
464
465     call  mix_SIA_L1 
466
467  end if spinup_4_vitdyn
468
469end subroutine step_thermomeca
Note: See TracBrowser for help on using the repository browser.