source: trunk/SOURCES/steps_time_loop.f90 @ 23

Last change on this file since 23 was 19, checked in by dumas, 9 years ago

climat-forcage-insolation_mod.f90 (Methode JB) validee avec fichier Ning, testsort_time avec passage de la variable time en double precision pour les sorties

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