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

Last change on this file since 40 was 25, checked in by dumas, 9 years ago

Suppression du module inutilise prop-thermiques_mod.f90

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