source: trunk/SOURCES/Old-sources/step.F90 @ 4

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

initial import GRISLI trunk

File size: 13.1 KB
Line 
1!> \file step.F90
2!! TOOOOOOOOO DOOOOOOOOOOOOO
3!! Step_grisli contient
4!<
5
6
7subroutine step_grisli()
8
9  use module3d_phy
10  use module_choix ! module de choix du type de run
11  !  module_choix donne acces a tous les modules
12  !  de declaration des packages
13  use interface_icetempmod
14  use sorties_ncdf_grisli
15  use flottab_mod
16  use diagno_mod
17  use resolmeca_SIA_L1
18  use track_debug 
19
20  use geom_type
21  use temperature_type
22  use ice_flow_type
23  use mask_flgz_type
24  use deformation_type
25  use autre_pr_temp_type
26  use step_temp_type
27 
28 
29implicit none
30
31  integer :: m
32  logical :: shelf_vitbil = .true.   ! si vrai on prend les vitesses de bilan pour le shelf
33
34!!$if (time.lt.timemax-100) then
35!!$   print*,'le modele remonte le temps de',timemax,'a',time
36!!$   stop
37!!$endif
38  timemax=time
39
40  if (itracebug.eq.1) write(num_tracebug,*) 'entree dans step_grisli', &
41       '   nt=',nt,'iter_beta=',iter_beta
42  !--------------------------------------------------------------------
43  !     lecture dans un petit fichier au cas ou les parametres du run 
44  !     seraient changes                                               
45  !--------------------------------------------------------------------
46
47  call limit_file(nt,real(time),dt,tend,dtsortie,dtcpt,testdiag,dtt,runname)
48
49  if (time.lt.tgrounded) then
50     marine=.false.
51  else
52     marine=.true.
53  endif
54
55  !if ((iter_beta.eq.0).or.(nt.le.1)) then
56  call masque()
57  call flottab()
58  if (iter_beta.gt.0)   call init_dragging
59  !end if
60
61
62  !-----------------------------------------------------------------------!
63  !              debut bloc appels tous les dtt                           !
64  !-------------   -------------------------------------------------------!
65
66
67  isync_1: if (ISYNCHRO.eq.1) then          !debut bloc appels tous les dtt   
68
69
70
71     spinup_1_veloc: if (ispinup.ne.2.or.nt.lt.5) then   
72
73     !---------------------------------------------------------------------------------
74     ! spinup_1 : calcul des vitesses dynamiques :
75     !------------
76     !            - est fait dans tous les cas pour les premieres boucles temporelles
77     !            - est fait la plupart des types de spinup sauf ispinup=2
78     !            - ispinup=1 ->  modif: maintenant on fait le calcul de veloc
79     !---------------------------------------------------------------------------------
80
81        if (itracebug.eq.1)  call tracebug(' Dans spinup_1_veloc : calculer velocities')
82
83        call forclim()                                  !
84        call ablation
85
86        call Neffect()                                             
87        call velocities()                                           
88     endif spinup_1_veloc
89
90     if (nt.gt.2) then
91        ! Hassine ca va etre modifier bien sur apres: l appel sera beaucoup moin long que ca
92
93        call Allocate_geom(geom_g,nx,ny)
94        call Init_geom(Geom_g)               
95         
96        call Allocate_temperature(temperature_g,nx,ny,nz,nz+nzm)
97        call Init_Temperature(Temperature_g)
98
99        call Allocate_Ice_flow(Ice_flow_g,nx,ny,nz)
100        call Init_Ice_flow(Ice_flow_g)
101       
102        call Allocate_Mask_flgz(Mask_flgz_g,nx,ny)
103        call Init_Mask_flgz(Mask_flgz_g)
104       
105        call Allocate_Deformation(Deformation_g,nx,ny,nz,n1poly,n2poly)
106        call Init_Deformation_h(Deformation_g)
107       
108        call Allocate_autre_pr_temp(Autre_pr_temp_g,nx,ny)
109        call Init_autre_pr_temp(Autre_pr_temp_g)
110     
111        call Init_step_temp(step_temp_g)
112       
113        call icetemp(geom_g,temperature_g, &
114             Ice_flow_g,Mask_flgz_g,&
115             deformation_g,Autre_pr_temp_g,&
116             step_temp_g,num_tracebug)
117
118        T= Temperature_g%Temperature 
119        Tpmp= Temperature_g%Temp_melting
120        Tbdot= Autre_pr_temp_g%Tbdot
121        Bmelt= Autre_pr_temp_g%BMELT
122        Ibase=Autre_pr_temp_g%IBASE
123        Phid=Autre_pr_temp_g%Phid
124
125        !      write(6,*) time,'call icetemp'
126        ! subroutines pour le calcul de la fusion basale
127        call bmeltshelf
128        call bmelt_grounded
129
130
131     endif
132
133    spinup_2_flowlaw: if ((ispinup.ne.2).or.(nt.lt.5)) then
134
135    !---------------------------------------------------------------------------------
136    ! spinup_2_flowlaw
137    !----------------
138    ! si ispinup =   2 , on ne fait pas le calcul de flowlaw, ni de diffusiv
139    !                    puisque les vitesses de bilan sont imposees (initialisee nt<5)
140    !                   
141    ! Flowlaw et diffusiv : servent soit directement (ispinup = 0,1) soit pour le
142    !                       rescaling (ispinup = 3)
143    !---------------------------------------------------------------------------------
144
145        if (itracebug.eq.1)  call tracebug(' Dans spinup_2_flowlaw')
146
147        if (ispinup.eq.3) then
148           call diffusiv                  ! pour vérifier les champs
149           debug_3D(:,:,71)=ux(:,:,1)
150           debug_3D(:,:,72)=uy(:,:,1)
151           debug_3D(:,:,13)=ux(:,:,nz)
152           debug_3D(:,:,14)=uy(:,:,nz)
153        end if
154
155        call flow_general
156        do iglen=n1poly,n2poly
157           call flowlaw(iglen)
158        end do
159
160
161!        debug_3D(:,:,65)=s2a_mx(:,:,1,1)
162!        debug_3D(:,:,66)=s2a_mx(:,:,1,2)
163!        debug_3D(:,:,67)=s2a_my(:,:,1,1)
164!        debug_3D(:,:,68)=s2a_my(:,:,1,2)
165!        debug_3D(:,:,69)=ddx(:,:,1)
166!        debug_3D(:,:,69)=ddx(:,:,2)
167
168        if (ispinup.eq.3) then
169           call calc_coef_vitbil                    ! recalibre sa_mx et s2a_mx et ddbx
170           !         call limit_coef_vitbil         ! garde la loi de deformation (debug)
171        end if
172
173     end if spinup_2_flowlaw
174
175  endif isync_1
176
177
178  spinup_3_bed: if (ispinup.eq.0.or.nt.lt.5) then
179    !---------------------------------------------------------------------------------
180    ! spinup_3_bed
181    !----------------
182    ! Ne passe dans ce bloc que quand on veut l'isostasie, donc si variations epaisseur
183    ! run standard  (ispinup=0) et pour les initialisations                 
184    !---------------------------------------------------------------------------------
185
186     if (itracebug.eq.1)  call tracebug(' Dans spinup_3_bed')
187
188     if ((nbed.eq.1).and.nt.ne.1.and.isynchro.eq.1) then           
189        call bedrock()                                            !  bedrock adjustment
190     endif
191  end if spinup_3_bed
192
193 spinup_4_vitdyn: if  (ispinup.ge.1.or.nt.lt.5.) then  ATTENTION A MODIFIER ipsinup=3  shelf_vitbil
194
195    !---------------------------------------------------------------------------------
196    ! spinup_4_vitdyn
197    !----------------
198    ! Ne passe dans ce bloc que quand on veut calculer les vitesses dynamiques :
199    ! - initialisations (nt.lt.5)
200    ! - run standard  (ispinup=0) (y compris iter_beta.ne.0)
201    ! - spinup temperature avec vitesses dynamiques  (ispinup = 1)               
202    !---------------------------------------------------------------------------------
203
204     !     (*********** DIFFUSIVITIES - VITESSES SHELF **********************)
205
206
207     call Neffect()
208
209     call diffusiv()
210
211
212iterbeta:     if (iter_beta.ne.0) then          ! pour les iterations sur beta
213
214        if (.not.shelf_vitbil) then    !  bloc si les vitesses shelves viennent de diagnoshelf
215
216                                       !  Vcol declare dans le module spinup_mod
217                                       !  ou dans le dragging si no_spinup
218
219           where (flotmx(:,:))
220              uxbar(:,:)=uxflgz(:,:)
221           elsewhere
222              uxbar(:,:)=Vcol_x(:,:)
223           end where
224
225           where (flotmy(:,:))
226              uybar(:,:)=uyflgz(:,:)
227           elsewhere
228              uybar(:,:)=Vcol_y(:,:)
229           end where
230        else
231           uxbar(:,:)=Vcol_x(:,:)
232           uybar(:,:)=Vcol_y(:,:)
233        end if
234     end if iterbeta
235
236
237     call strain_rate()
238 
239    if (iter_beta.ne.0) then
240        uxbar(:,:)=uxflgz(:,:)+uxdef(:,:)
241        uybar(:,:)=uyflgz(:,:)+uydef(:,:)
242     end if
243
244
245
246     !-----------------------------------------------------------------------
247     !              debut bloc appels tous les dtt                         
248     !-------------   -------------------------------------------------------
249
250     if ((shelfy).and.(marine).and.(isynchro.eq.1)) then               
251
252        call lake_to_slv
253        !write(6,*) 'premier diagno dans step'
254        !call flottab
255
256        ! les 3 lignes ci-dessous pour avoir les champs avant diagno (en cas de pb colineaire)
257        !call init_sortie_ncdf
258        !call sortie_ncdf_cat
259
260        call diagnoshelf
261
262        !call sortie_ncdf_cat
263        ! STOP                 
264
265
266        !         faire 10 tours pour equilibrer                               
267        if ((icompteur.eq.0).and.(nt.lt.10)) then                       
268           do m=1,1                                   
269              !         write(6,*) 'dans boucle m'
270              call diagnoshelf()                                       
271           end do
272        end if
273
274     endif !
275
276     call  mix_SIA_L1
277
278     debug_3D(:,:,78)= abs(uxbar(:,:))-abs(Vcol_x(:,:))
279     debug_3D(:,:,79)= abs(uybar(:,:))-abs(Vcol_y(:,:))
280
281     debug_3D(:,:,80)= abs(uxdef(:,:))-abs(Vcol_x(:,:))
282     debug_3D(:,:,81)= abs(uydef(:,:))-abs(Vcol_y(:,:))
283
284
285  end if spinup_3_bed
286
287end subroutine step_grisli
288
289!--------------------------------------------------------------------------
290! step_grisli1 est a l'interieur de la boucle temps.
291! il commence par icethick, qui détermine dt, isynchro
292
293! ensuite contient les appels aux diverses sorties
294! puis un appel a step_grisli
295
296
297subroutine step_grisli1               
298
299  use module3d_phy
300
301  use geom_type
302  use temperature_type
303  use ice_flow_type
304  use mask_flgz_type
305  use deformation_type
306  use autre_pr_temp_type
307
308  use module_choix ! module de choix du type de run
309  !  module_choix donne acces a tous les modules
310  !  de declaration des packages
311  use sorties_ncdf_grisli
312  use flottab_mod
313  use interface_icetempmod
314  use diagno_mod 
315  use track_debug 
316
317  if (itracebug.eq.1) write(num_tracebug,*) ' Entree dans step_grisli1',   &
318       '   nt=',nt,'iter_beta=',iter_beta
319
320  nospinup_4: if ((ispinup.eq.0.or.nt.lt.5).and.(iter_beta.eq.0)) then
321     if (itracebug.eq.1)  call tracebug(' Dans nospinup_4 : appelle icethick')
322
323     call icethick3()
324
325     call calving
326
327     if (itracebug.eq.1)  call tracebug('apres calving')
328
329  else if ((ispinup.eq.1).or.(ispinup.eq.3).or.(iter_beta.gt.0)) then
330     dt=dtmax
331     hdot(:,:)=0.
332
333     !      if (nt.eq.6) time=0.
334     !      time=nint(dtmax*nint(time/dtmax))+dtmax  semble poser des problemes
335     time=nt*dtmax
336
337     if (itracebug.eq.1) write(num_tracebug,*) ' Dans nospinup_4  no icethick ',ispinup, &
338          '    time=',time
339  else if (ispinup.eq.2) then
340     call force_balance_vel
341     call icethick3()
342     call calving
343
344  end if nospinup_4
345
346  CALL write_trace_940('ICE SHEET SURFACE')
347
348  where(.not.flot(:,:))
349     S(:,:)=Bsoc(:,:)+H(:,:)
350     B(:,:)=Bsoc(:,:)
351  end where
352
353  ! calcul de Hmx et Hmy -> shift=-1, dim=1 -> H(i-1,j)
354
355  hmx(:,:)=0.5*(H(:,:)+eoshift(H(:,:),shift=-1,boundary=0.,dim=1))
356  hmy(:,:)=0.5*(H(:,:)+eoshift(H(:,:),shift=-1,boundary=0.,dim=2))
357  hmx(1,:)=0.
358  hmy(:,1)=0.
359
360  ! mise a 0 des vitesses quand epaisseur nulle
361  where (hmx(:,:).le.1.)
362     uxbar(:,:)=0.
363  endwhere
364  where (hmy(:,:).le.1.)
365     uybar(:,:)=0.
366  endwhere
367
368  !     ________________________________  fin bloc appel tous les tend
369
370  !     (****** ICE SHEET STATISTICS and SHORT OUTPUT *******)
371  !     if (mod(nt,nint(NDISP/DT)).eq.0).or.(nt.eq.1).or.(nt.eq.2).or.
372  !    &  (nt.eq.3).or.(nt.eq.4).or.(nt.eq.5).or.(nt.eq.10).or.
373
374  !     if ((mod(abs(TIME),NDISP*1.).lt.dtmin).or.(NT.le.5).or.(NT.eq.10) &
375  !            .or.(nt.eq.15).or.(nt.eq.20).or.(nt.eq.30).or. &
376  !            (abs(TIME-TEND).lt.dtmin)) then
377
378  !        call sortie_ncdf_cat
379
380
381  if ((mod(abs(TIME),NDISP*1.).lt.dtmin).or.(isynchro.eq.1)) then
382
383     call shortoutput()
384
385     if ((geoplace.eq.'anteis1') &
386          .or.(geoplace.eq.'ant20km')) then
387
388        call ts_out()
389     endif
390     !        call sealevel_out
391
392  endif
393
394
395  ! Sorties horizontales.
396  !--------------------------------------------------------------------------------
397  !  Premiere partie des sorties horizontales
398
399  !     if (mod(abs(dble(TIME)),dble(DTSORTIE)).lt.dble(dtmin)) then
400  !
401
402  if ((mod(abs(TIME),dtt).lt.dtmin).or.(isynchro.eq.1)) then
403     isynchro=1
404     !     if (isynchr o.eq.1) then
405     call testsort_time(real(time))
406     if (iglob_hz.eq.1) then
407        call sortie_hz_multi   ! pour les variables déclarées dans 3D
408        !              call hz_output(real(time))
409     endif
410
411  else
412     iglob_hz=0
413
414  endif
415
416
417  !   *** Sortie des profils tous les dtprofile ***
418  if ((NT.eq.1) &    !   .or.(NT.eq.2)
419       .or.(mod(abs(dble(TIME)),dble(DTPROFILE)).lt.dble(dtmin)) &
420       .or.(abs(TIME-TEND).lt.dtmin)) then
421
422     !        call sortieprofile()
423
424  endif
425
426  !     *** SORTIE COMPTEUR TOUS LES DTCPT ANS ***
427
428  if ((mod(abs(dble(TIME)),dble(DTCPT)).lt.dble(dtmin)).OR.   &
429       (ABS(TIME+297000).LT.dtmin).OR.   &
430       (ABS(TIME+126000).LT.dtmin).OR.   &
431       (ABS(TIME+21000).LT.dtmin).OR.    &
432       (ABS(TIME+16000).LT.dtmin))    THEN
433
434     call sortiecptr()
435
436  endif
437
438
439  !  ************* sorties netcdf ****************
440
441  call testsort_time_ncdf(real(time))
442  if (iglob_ncdf .EQ. 1) call sortie_ncdf_cat
443
444  !   call sortie_ncdf_cat
445
446  call track_change_T
447
448  call step_grisli()
449  !
450  !
451  !
452end subroutine step_grisli1
Note: See TracBrowser for help on using the repository browser.