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

Last change on this file since 52 was 52, checked in by dumas, 8 years ago

Merge branche iLOVECLIM sur rev 51

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