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