source: trunk/SOURCES/steps_time_loop.f90 @ 170

Last change on this file since 170 was 165, checked in by dumas, 7 years ago

Spinup with iterbeta is now compatible with steps_time_loop.f90
File size: 12.8 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 bilan_flux_mod
26  !  use track_debug
27
28  implicit none
29  integer :: nt_init = 0      ! number of loops for initialisation of ice thickness
30  integer :: iter_visco       ! number of iterations for ssa viscosity
31
32  if (itracebug.eq.1)  call tracebug('Entree dans step_time_loop ')
33  if (itracebug.eq.1) write(num_tracebug,*) 'nt = ',nt
34
35
36  ! ice thickness evolution
37  !=============================
38
39
40  spinup_icethick: if (ispinup.eq.0.or.ispinup.eq.2.or.nt.lt.nt_init) 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)
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     call flottab
59     call determin_tache
60     call calving
61     call ablation_bord
62     call bilan_eau
63     call bilan_flux_output     
64     call flottab
65     if (isynchro.eq.1.or.nt.eq.1) call determin_marais
66     if (itracebug.eq.1)  call tracebug('apres marais')
67
68     ! new calculation of ice surface elevation
69     where(.not.flot(:,:))
70        S(:,:)=Bsoc(:,:)+H(:,:)
71        S(:,:)=max(S(:,:),sealevel)
72        B(:,:)=Bsoc(:,:)
73     end where
74
75     ! calcul de Hmx et Hmy -> shift=-1, dim=1 -> H(i-1,j)
76
77     hmx(:,:)=0.5*(H(:,:)+eoshift(H(:,:),shift=-1,boundary=0.,dim=1))
78     hmy(:,:)=0.5*(H(:,:)+eoshift(H(:,:),shift=-1,boundary=0.,dim=2))
79     hmx(1,:)=0.
80     hmy(:,1)=0.
81
82     ! mise a 0 des vitesses quand epaisseur nulle
83     where (hmx(:,:).le.1.)
84        uxbar(:,:)=0.
85     endwhere
86     where (hmy(:,:).le.1.)
87        uybar(:,:)=0.
88     endwhere
89
90  else   ! time increment when no icethick : ispinup= 1 or 3    ! loop spinup_icethick
91
92     hdot(:,:)=0.
93     time=nt*dtmax
94
95     if (itracebug.eq.1) write(num_tracebug,*) ' Dans spinup_icethick:  no icethick ',ispinup, &
96          '    time=',time
97
98
99  end if spinup_icethick
100  ! end of ice thickness evolution
101  !====================================================================
102
103
104
105  ! outputs
106  !====================================================================
107
108  !   vertical plan snapshots (profiles)
109  !------------------------------------------------------------------
110
111  !if ((NT.eq.1) &    !   .or.(NT.eq.2)
112  !     .or.(mod(abs(dble(TIME)),dble(DTPROFILE)).lt.dble(dtmin)) &
113  !     .or.(abs(TIME-TEND).lt.dtmin)) then
114  !  call sortieprofile()
115  !endif
116
117  !  sorties netcdf Hassine (aout 2010)  (2D and 3D)
118  !------------------------------------------------------------------
119
120  call testsort_time_ncdf(time)
121  if (iglob_ncdf .EQ. 1) call sortie_ncdf_cat
122  if ((isynchro.eq.1).or.(nt.eq.1)) then
123    call shortoutput
124    diff_H = 0.
125    Bm_dtt(:,:) = 0.
126    bmelt_dtt(:,:) = 0.
127    calv_dtt(:,:)=0.
128    ablbord_dtt(:,:)=0.
129    diff_H_2D(:,:)=0.
130    grline_dtt(:,:)=0.
131  endif
132
133  !   sortie compteur tous les dtcpt ans
134  !------------------------------------------------------------------
135  !iout == 1 sortie cptr
136  !iout == 2 sortie nc
137  if (itracebug.eq.1)  call tracebug('dans steps_time_loop avant call out_recovery ')
138
139  call out_recovery(iout)
140 
141  ! end of outputs
142
143  !====================================================================
144  !
145  call step_thermomeca()
146  !
147  ! contient tout le reste de la dynamique et du forçage
148  !====================================================================
149
150end subroutine step_time_loop
151
152
153
154
155
156
157!---------------------------------------------------------------------
158! The subroutine  step_thermomeca() calls the following processes
159!
160! - climatic forcing and its impact on surface mass balance
161! - dynamic velocities (SIA)
162! - ice temperature
163
164subroutine step_thermomeca()
165
166
167  use module3d_phy
168  use module_choix ! module de choix du type de run
169  !  module_choix donne acces a tous les modules
170  !  de declaration des packages
171  use icetempmod
172  use sorties_ncdf_grisli
173  use flottab_mod
174  use diagno_mod
175  use resolmeca_SIA_L1
176!  use track_debug
177
178
179  implicit none
180
181  integer :: m
182  logical :: shelf_vitbil = .true.   ! si vrai on prend les vitesses de bilan pour le shelf
183  integer :: nt_init_tm = 0          ! number of loops for initialisation of thermo mechanical
184  integer :: iter_visco              ! number of iterations for ssa viscosity
185
186  if (ispinup.le.1) shelf_vitbil = .false.         ! general case, ice shelves velocities are computed by diagnoshelf
187
188  timemax=time
189  if (itracebug.eq.1) write(num_tracebug,*) 'entree dans step_thermomeca', &
190       '   nt=',nt
191
192  !--------------------------------------------------------------------
193  !     lecture dans un petit fichier "runname.limit" au cas ou les
194  !     parametres du run seraient changes 
195  !                                                   
196  !--------------------------------------------------------------------
197
198!  call limit_file(nt,real(time),dt,tend,dtsortie,dtcpt,testdiag,dtt,runname)
199
200  !if (time.lt.tgrounded) then
201  !   marine=.false.
202  !else
203  marine=.true.
204  !endif
205
206  ! update the regions (floating, ice ...)
207  !-------------------------------------------
208  call masque()
209
210
211  !-----------------------------------------------------------------------!
212  !   debut d'un bloc appels tous les dtt (pas de temps long- asynchrone) !
213  !-------------   -------------------------------------------------------!
214
215
216
217  isync_1: if (ISYNCHRO.eq.1) then          ! debut bloc appels tous les dtt   
218
219     ! climatic forcing
220     !=====================
221     call forclim                             
222     call ablation
223
224
225
226     ! SIA dynamic velocities  (deformation + loi de glissement )
227     !========================
228
229     ! remark : velocity is 3D
230     !          to update the velocity for ice temperature calculation
231     !            - surface mass balance for velocity (because ablation was called)
232     !            - impact of surface slope
233
234
235     spinup_1_veloc: if ((ispinup.ne.2).or.(nt.lt.nt_init_tm)) then        ! cas general
236
237        !---------------------------------------------------------------------------------
238        ! spinup_1 : calcul des vitesses dynamiques :
239        !------------
240        !            - est fait dans tous les cas pour les premieres boucles temporelles
241        !            - est fait pour ispinup=3 pour calibration forme de vitesse
242        !            - n'est pas fait pour ispinup=2 (test balance vel sur icethick)
243        !
244        !---------------------------------------------------------------------------------
245
246        if (itracebug.eq.1)  call tracebug(' Dans spinup_1_veloc : calculer SIA_velocities')
247
248        call Neffect()                                             
249        call SIA_velocities()
250
251     else    ! ispinup = 2
252
253        call force_balance_vel   ! fait l'ensemble du champ (flgz ou pas)
254
255     endif spinup_1_veloc
256
257
258     ! ice temperature
259     !===================================================================
260
261     ! update values in the structures Geom_g, Temperature_g, ...
262
263! afq     calc_temp:     if ((nt.gt.2).and.(geoplace(1:5).ne.'mism3')) then
264     calc_temp:     if (geoplace(1:5).ne.'mism3') then
265        if (itracebug.eq.1)  call tracebug('avant appel icetemp')
266        call icetemp
267        ! subroutines pour le calcul de la fusion basale
268        call bmeltshelf
269        call bmelt_grounded
270
271     endif calc_temp
272
273
274     ! flowlaw
275     !===================================================================
276
277     spinup_2_flowlaw: if ((ispinup.ne.2).or.(nt.lt.nt_init_tm)) then    ! cas general
278
279        !---------------------------------------------------------------------------------
280        ! spinup_2_flowlaw
281        !----------------
282        ! si ispinup =   2 , on ne fait pas le calcul de flowlaw, ni de diffusiv
283        !                    puisque les vitesses de bilan sont imposees
284        !                                             (initialisee nt < nt_init_tm)
285        !                   
286        ! Flowlaw et diffusiv : servent soit directement (ispinup = 0,1) soit pour le
287        !                       rescaling (ispinup = 3)
288        !---------------------------------------------------------------------------------
289
290        if (itracebug.eq.1)  call tracebug(' Dans spinup_2_flowlaw')
291
292           call flow_general
293           do iglen=n1poly,n2poly
294              call flowlaw(iglen)
295           end do
296
297
298     end if spinup_2_flowlaw
299
300     if (ispinup.eq.3) then
301        call calc_coef_vitbil                    ! recalibre sa_mx et s2a_mx et ddbx pour toute la zone posee
302     end if                                      ! Dans le shelf flgz = Ubar = Vbil et udef=0
303
304     call tracer      ! aurel neem (attention, position a verifier)   
305
306  endif isync_1
307
308
309  ! fin du bloc tous les dtt isync_1
310  !------------------------------------
311
312
313  spinup_4_vitdyn: if  ((ispinup.le.1).or.(nt.lt.nt_init_tm)) then 
314
315     !---------------------------------------------------------------------------------
316     ! spinup_4_vitdyn
317     !----------------
318     ! Passe dans ce bloc que quand on veut calculer les vitesses dynamiques :
319     ! - initialisations (nt.lt.nt_init_tm)
320     ! - run standard  (ispinup=0)
321     ! - spinup temperature avec vitesses dynamiques  (ispinup = 1)               
322     !---------------------------------------------------------------------------------
323
324
325
326     call Neffect()
327
328     call diffusiv()
329
330     ! strain rate calculations
331     !===========================
332
333     call strain_rate()
334
335     ! ice shelves and ice streams velocities : SSA
336     !==============================================
337
338     !-----------------------------------------------------------------------
339     !              debut d'un deuxieme bloc appels tous les dtt                         
340     !-------------   -------------------------------------------------------
341
342     isync_2:    if ((shelfy).and.(marine).and.(isynchro.eq.1)) then               
343
344        !les 3 lignes ci-dessous pour avoir les champs avant diagno (en cas de pb colineaire)
345        !call init_sortie_ncdf
346        !call sortie_ncdf_cat
347
348        if (icompteur.eq.1) then    ! reprise de tout le vecteur d'etat   
349           iter_visco = 1
350
351        else if ((icompteur.ne.1).and.(nt.le.1)) then  ! initialisation avec reprise partielle
352                                                       ! ou nulle du vecteur d' etat
353          iter_visco = 10
354        else
355          iter_visco = 2
356        end if
357 
358           test_iter_diag=0.
359           do m=1,iter_visco             
360              call diagnoshelf
361              call mix_SIA_L1 
362              call strain_rate
363              if (test_iter_diag.lt.1.e-2) exit
364           end do
365           
366           ! iterations pour inversion du beta en fonction des vitesses de bilan
367           if ((inv_beta.eq.1).and.(time.ge.time_iter).and.(time.le.time_iter_end)) then  ! iterations pendant 20 ans
368              do m=1,nb_iter_vitbil
369                 call beta_iter_vitbil(m)
370                 call diagnoshelf
371                 call mix_SIA_L1
372                 call strain_rate
373              end do
374           end if
375           ! fin des iterations beta_iter_vitbil
376
377        endif isync_2
378
379     call  mix_SIA_L1 
380
381  end if spinup_4_vitdyn
382 
383 
384  ! isostasy
385  !=================
386
387  spinup_3_bed: if ((ISYNCHRO.eq.1).and.(ispinup.eq.0.or.nt.lt.nt_init_tm)) then
388  !---------------------------------------------------------------------------------
389  ! spinup_3_bed
390  !----------------
391  ! Ne passe dans ce bloc que quand on veut l'isostasie, donc si variations epaisseur
392  ! run standard  (ispinup=0) et pour les initialisations                 
393  !---------------------------------------------------------------------------------
394
395  if (itracebug.eq.1)  call tracebug(' Dans spinup_3_bed')
396
397  if ((nbed.eq.1).and.nt.ne.1.and.isynchro.eq.1.and.(mod(abs(TIME),50.).lt.dtmin)) then
398    call bedrock                 !  bedrock adjustment
399  endif
400
401  call lake_to_slv               ! pour traiter dynamiquement les shelves sur lac
402
403  end if spinup_3_bed
404 
405
406end subroutine step_thermomeca
Note: See TracBrowser for help on using the repository browser.