source: branches/iLoveclim/SOURCES/branche-Cat-spinup-dec2011/steps_time_loop.f90 @ 77

Last change on this file since 77 was 4, checked in by dumas, 10 years ago

initial import GRISLI trunk

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