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

Last change on this file since 100 was 91, checked in by aquiquet, 8 years ago

GRISLI coupled using basal melting rates coming from CLIO. Might need further testing for LGM conditions.

File size: 14.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 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     ! oceanic forcing
249     !=====================
250     call bmeltshelf  ! afq -- icetemp is not called before the 4th dtt, so a first call to bmelt here
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        ! afq -- 28-10-2016: there is a reset of bmelt in temp_col in the
294        ! current version of the code, so shelf and grounded melt have to
295        ! be called here. TODO: remove reset in icetemp routines??
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         if ((nbed.eq.1).and.nt.ne.1.and.isynchro.eq.1.and.(mod(abs(TIME),50.).lt.dtmin)) then
349           call bedrock                                              !  bedrock adjustment
350        endif
351
352        call lake_to_slv                                             ! pour traiter dynamiquement les shelves sur lac
353
354     end if spinup_3_bed
355
356     call tracer      ! aurel neem (attention, position a verifier)   
357
358  endif isync_1
359
360
361  ! fin du bloc tous les dtt isync_1
362  !------------------------------------
363
364
365  spinup_4_vitdyn: if  ((ispinup.le.1).or.(nt.lt.nt_init_tm)) then 
366
367     !---------------------------------------------------------------------------------
368     ! spinup_4_vitdyn
369     !----------------
370     ! Passe dans ce bloc que quand on veut calculer les vitesses dynamiques :
371     ! - initialisations (nt.lt.nt_init_tm)
372     ! - run standard  (ispinup=0) (y compris iter_beta.ne.0)
373     ! - spinup temperature avec vitesses dynamiques  (ispinup = 1)               
374     !---------------------------------------------------------------------------------
375
376
377
378     call Neffect()
379
380     call diffusiv()
381
382
383     iterbeta:     if (iter_beta.ne.0) then          ! only for iterations on beta
384
385        if (.not.shelf_vitbil) then    !  bloc si les vitesses shelves viennent de diagnoshelf
386
387           !  Vcol declare dans le module spinup_mod
388           !  ou dans le dragging si no_spinup
389
390           where (flotmx(:,:))
391              uxbar(:,:)=uxflgz(:,:)
392           elsewhere
393              uxbar(:,:)=Vcol_x(:,:)
394           end where
395
396           where (flotmy(:,:))
397              uybar(:,:)=uyflgz(:,:)
398           elsewhere
399              uybar(:,:)=Vcol_y(:,:)
400           end where
401        else
402           uxbar(:,:)=Vcol_x(:,:)
403           uybar(:,:)=Vcol_y(:,:)
404        end if
405     end if iterbeta
406
407
408     ! strain rate calculations
409     !===========================
410
411     call strain_rate()
412
413
414     if (iter_beta.ne.0) then               ! on n'a utilise les vitesses de bilan que pour strain_rate "Dirichlet ?"
415        uxbar(:,:)=uxflgz(:,:)+uxdef(:,:)
416        uybar(:,:)=uyflgz(:,:)+uydef(:,:)
417     end if
418
419
420
421     ! ice shelves and ice streams velocities : SSA
422     !==============================================
423
424     !-----------------------------------------------------------------------
425     !              debut d'un deuxieme bloc appels tous les dtt                         
426     !-------------   -------------------------------------------------------
427
428     isync_2:    if ((shelfy).and.(marine).and.(isynchro.eq.1)) then               
429
430        !les 3 lignes ci-dessous pour avoir les champs avant diagno (en cas de pb colineaire)
431        !call init_sortie_ncdf
432        !call sortie_ncdf_cat
433
434        if (icompteur.eq.1) then    ! reprise de tout le vecteur d'etat   
435           iter_visco = 1
436
437        else if ((icompteur.ne.1).and.(nt.le.1)) then  ! initialisation avec reprise partielle
438                                                       ! ou nulle du vecteur d' etat
439          iter_visco =  10
440        else
441           iter_visco = 2
442        end if
443 
444        iter_visco= 10                                 ! warning, test sur l'impact iteration dragging et visco
445                                                       ! Cat 18 janv 2013
446
447           do m=1,iter_visco             
448              call diagnoshelf
449              call  mix_SIA_L1 
450              call strain_rate()
451              call flottab()                           ! pour iteration dragging : flottab appelle dragging
452!              write(444,*) time,test_iter_diag,m
453!              if (test_iter_diag.lt.1.e-4) exit
454!              write(555,*) time,test_iter_diag,m
455!               write(222,*) time,test_iter_diag,m
456              if (test_iter_diag.lt.1.e-2) exit
457!              write(333,*) time,test_iter_diag,m
458!              if (test_iter_diag.lt.1.e-3) exit
459
460           end do
461        endif isync_2
462
463
464     call  mix_SIA_L1 
465
466  end if spinup_4_vitdyn
467
468end subroutine step_thermomeca
Note: See TracBrowser for help on using the repository browser.