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

Last change on this file since 151 was 151, checked in by aquiquet, 7 years ago

Grisli-iloveclim branch merged to trunk at revision 150 + bug correction for surface temperature in climat_coupl

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