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