source: trunk/SOURCES/Draggings_modules/beta_iter_vitbil_mod.f90 @ 70

Last change on this file since 70 was 70, checked in by dumas, 8 years ago

Bug fix in Qprod_icetemp (instability in ice temperature) : slowssa points heat = deformation SIA + sliding. Use fleuvemx and fleuvemy to identifie streams SSA and slow SSA points
File size: 16.4 KB
Line 
1!< \file beta_iter_vitbil_mod.f90
2!< Module pour affiner le beta a partir d'une carte vitbil  sur les noeuds centres.
3!< A appeler depuis steps_time_loop
4!< a utiliser avec beta_prescr
5
6!> \ beta_iter_vitbil
7!! \author ...
8!! \date ...
9!! @note Used module
10!! @note   - use module3D_phy
11!<
12
13module beta_iter_vitbil_mod
14
15  use module3d_phy
16  use interface_input
17  use netcdf
18  use io_netcdf_grisli
19
20  implicit none
21
22
23  real*8, dimension(:,:),   pointer  :: tab => null() !< tableau 2d real ecrit dans le fichier pour lect ncdf
24
25  real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite
26  real, dimension(nx,ny) :: Vcol_y           !<
27  real, dimension(nx,ny) :: Vsl_x            !<
28  real, dimension(nx,ny) :: Vsl_y            !<
29
30
31  real, dimension(nx,ny) :: drag_centre      !< beta on major node (average)
32  real :: beta_limgz                         !< when beta gt beta_limgz -> not gzmx
33  real :: beta_min                           !< minimum value of beta
34  real :: beta_mult                          !< coefficient of beta field
35  integer  :: ill,jll,nmoy
36
37  real :: maxi                               !< calcul correction dS a la grounding line
38  real :: mini
39
40  logical :: corr_def = .False.               !< for deformation correction (compatibility)
41
42
43  real, dimension(nx,ny) :: Umag_vitbil             !< vitesse de bilan
44  real, dimension(nx,ny) :: Uslid_vitbil            !< vitesse de bilan partie glissement
45  real, dimension(nx,ny) :: Umag_direct             !< vitesse moyenne centree calculee par GRISLI
46  real, dimension(nx,ny) :: Uslmag_direct           !< vitesse glissement centree calculee par GRISLI
47  real, dimension(nx,ny) :: Udefmag_direct          !< vitesse deformation centree calculee par GRISLI
48  real, dimension(nx,ny) :: driving_stress          !< driving stress
49
50  character(len=100)     :: Umag_bil_file           !< fichier qui contient les donnees
51  real                   :: time_iter               !< temps de debut des iterations
52  integer                :: nb_iter_vitbil          !< nombre d'iteratiosn avant de changer de pas de temps
53  real                   :: coef_iter_vitbil        !< coefficient <= 1 (rapport vitesses)
54  real                   :: J_Umag                  !< estimateur ecart entre vitesses
55  real                   :: J_Udef                  !< estimateur vitesse regions sans glissement
56  integer                :: nb_umag                 !< nombre de points calcul J_Umag
57  integer                :: nb_Udef                 !< nombre de points calcul J_Umag
58
59
60contains
61
62  !-----------------------------------------------------------------------------------------
63  !> SUBROUTINE: init_dragging
64  !! Cette routine fait l'initialisation du dragging.
65  !<
66  subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
67    ! nouvelle version : lit les fichiers nc
68    implicit none
69    character(len=100) :: beta_c_file               !  beta on centered grid
70    character(len=100) :: betamx_file               !  beta on staggered grid mx
71    character(len=100) :: betamy_file               !  beta on staggered grid mx
72
73    namelist/beta_prescr/beta_c_file,beta_limgz,beta_min,beta_mult
74
75    if (itracebug.eq.1)  call tracebug(' Prescr_beta       subroutine init_dragging')
76
77    iter_beta=0
78
79    ! lecture des parametres du run                   
80    !--------------------------------------------------------------------
81
82    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
83428 format(A)
84    read(num_param,beta_prescr)
85
86    write(num_rep_42,428) '!___________________________________________________________' 
87    write(num_rep_42,428) '!  read beta on centered grid'
88    write(num_rep_42,beta_prescr)                   
89    write(num_rep_42,428) '!  beta_file : nom des fichier qui contiennent les beta_c'
90    write(num_rep_42,428) '!  above beta_limgz, gzmx is false'
91    write(num_rep_42,428) '!  if grounded, beta_min minimum value '
92    write(num_rep_42,428) '!  beta_mult : coefficient just after reading' 
93    write(num_rep_42,*)
94    write(num_rep_42,428) '!___________________________________________________________' 
95
96
97    beta_c_file = trim(dirnameinp)//trim(beta_c_file)
98    betamx_file = trim(dirnameinp)//trim(betamx_file)
99    betamy_file = trim(dirnameinp)//trim(betamy_file)
100    betamax = beta_limgz
101    betamax_2d(:,:) = betamax
102
103!    call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'')
104    write(6,*)  beta_c_file
105
106    call Read_Ncdf_var('z',beta_c_file,tab)
107    drag_centre(:,:) = tab(:,:)
108
109
110    ! multiplication par le coefficient beta_mult
111
112    drag_centre(:,:) = drag_centre(:,:)*beta_mult
113
114    where (mk_init(:,:).eq.-2)                                             ! iles a probleme
115       drag_centre(:,:) = 2.* abs(beta_limgz)
116    end where
117
118    where(flot(:,:))
119       drag_centre(:,:) = 0.
120    end where
121
122
123    call  beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)         ! redistribute on staggered grid
124
125    if  (beta_limgz.gt.0.) then
126       beta_centre(:,:) =  drag_centre(:,:)
127       betamx(:,:)      =  drag_mx(:,:)
128       betamy(:,:)      =  drag_my(:,:)
129
130    else if (beta_limgz.lt.0.) then                              ! attention sans doute pour plastic
131
132       beta_centre(:,:) = - beta_limgz
133       betamx(:,:)      = - beta_limgz
134       betamy(:,:)      = - beta_limgz
135       drag_centre(:,:) = - beta_limgz 
136       drag_mx(:,:)     = - beta_limgz
137       drag_my(:,:)     = - beta_limgz 
138       beta_limgz = 1.
139    end if
140
141
142
143    call beta_min_value(beta_min)      !  valeur mini que peut avoir beta (en Pa an m-1)
144    ! on peut envisager une valeur ~ 10
145    ! rappel : beta doit être positif.
146
147
148
149    !    call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)                  ! centered distributed
150    !                                                                         ! on staggered
151    !
152
153
154    beta_centre(:,:) =  drag_centre(:,:)                                 ! surtout pour sorties
155
156
157
158    where (drag_mx(:,:).ge.beta_limgz)
159       mstream_mx(:,:) = 0
160       betamx(:,:)     = beta_limgz
161       drag_mx(:,:)    = beta_limgz
162    elsewhere
163       mstream_mx(:,:) = 1
164       betamx(:,:)     = drag_mx(:,:)
165    endwhere
166
167    where (drag_my(:,:).ge.beta_limgz)
168       mstream_my(:,:) = 0
169       betamy(:,:)     = beta_limgz
170       drag_my(:,:)    = beta_limgz
171    elsewhere
172       mstream_my(:,:) = 1
173       betamy(:,:)     = drag_my(:,:)
174    endwhere
175
176    if (itracebug.eq.1)  call tracebug(' Prescr_beta  apres lecture')
177    mstream_mx(:,:) = 0
178    mstream_my(:,:) = 0
179
180
181    call gzm_beta_prescr
182
183    call init_beta_iter_vitbil
184    return
185
186  end subroutine init_dragging
187  !________________________________________________________________________________
188
189!-----------------------------------------------------------------------------------------
190!< subroutine : init_beta_iter_vitbil
191!< Cette routine fait l'initialisation de la procedure beta_iter_vitbil
192
193subroutine init_beta_iter_vitbil
194
195
196  namelist/beta_iter_vitbil/time_iter,nb_iter_vitbil,coef_iter_vitbil,Umag_bil_file
197
198
199    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
200428 format(A)
201    read(num_param,beta_iter_vitbil)
202
203    write(num_rep_42,428) '!___________________________________________________________' 
204    write(num_rep_42,428) '!  module beta_iter_vitbil'
205    write(num_rep_42,beta_iter_vitbil)                   
206    write(num_rep_42,428) '! ' 
207    write(num_rep_42,428) '!  Umag_bil_file : fichier qui contient les donnees Umag_vitbil'
208    write(num_rep_42,428) '!  time_iter     : temps de debut des iterations '
209    write(num_rep_42,428) '!  nb_iter_vitbil: nombre diteratation a faire avant de changer de pas de temps ' 
210    write(num_rep_42,*)   '!  coef_iter_vitbil : coefficient pour le rapport vitesses  '
211    write(num_rep_42,428) '!___________________________________________________________' 
212
213    time_iter = time_iter + tbegin       ! time_iter est le temps de relaxation (de l'ordre de 5 ans)
214                                         ! ajouter tbegin pour ne pas dependre du temps de depart.
215
216
217    Umag_bil_file = trim(dirnameinp)//trim(Umag_bil_file)
218
219    write(6,*) Umag_bil_file
220!    call lect_input(1,'Umag_vitbil',1,Umag_vitbil,Umag_bil_file,'')
221    call Read_Ncdf_var('Umag_vitbil',Umag_bil_file,tab)
222    Umag_vitbil(:,:) = tab(:,:)
223   
224    debug_3D(:,:,58) = Umag_vitbil(:,:)
225
226  end subroutine init_beta_iter_vitbil
227
228
229  !________________________________________________________________________________
230
231  !> SUBROUTINE: dragging
232  !! Defini la localisation des streams et le frottement basal
233  !<
234  !-------------------------------------------------------------------------
235  subroutine dragging   ! defini la localisation des streams et le frottement basal
236
237    implicit none
238
239
240
241    if (itracebug.eq.1)  call tracebug(' beta_prescr         subroutine dragging')
242
243
244    betamx(:,:)=drag_mx(:,:)
245    betamy(:,:)=drag_my(:,:)
246
247
248    call  gzm_beta_prescr        ! determine gz et flgz et met a 0 le beta des points flottants
249
250
251    return
252  end subroutine dragging
253
254  !----------------------------------------------------------------------------------
255  !>SUBROUTINE: gzm_beta_prescr
256  !! Calcul de gzmx
257  !<
258
259  subroutine gzm_beta_prescr     ! attribue flgzmx et gzmx (et my)
260   
261    logical ::  test_Tx          ! test if there is a basal melting point in the surrounding
262    logical ::  test_Ty          ! test if there is a basal melting point in the surrounding
263    logical ::  test_beta_x      ! test if there is a low drag point in the surrounding
264    logical ::  test_beta_y      ! test if there is a low drag point in the surrounding
265
266    real    ::  beta_forc_fleuve ! below this value -> ice stream
267
268!    beta_forc_fleuve = 5.e3
269    beta_forc_fleuve = beta_limgz
270
271    ! calcul de gzmx
272
273
274    gzmx(:,:)=.true.
275    gzmy(:,:)=.true.
276    flgzmx(:,:)=.false.
277    flgzmy(:,:)=.false.
278
279
280    ! points flottants : flgzmx mais pas gzmx
281    do j=2,ny
282       do i=2,nx
283          if (flot(i,j).and.(flot(i-1,j))) then
284             flotmx(i,j)=.true.
285             flgzmx(i,j)=.true.
286             gzmx(i,j)=.false.
287             betamx(i,j)=0.
288          end if
289          if (flot(i,j).and.(flot(i,j-1))) then
290             flotmy(i,j)=.true.
291             flgzmy(i,j)=.true.
292             gzmy(i,j)=.false.
293             betamy(i,j)=0.
294          end if
295       end do
296    end do
297
298    if (itracebug.eq.1)  call tracebug(' apres criteres flot')
299
300    if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330) 
301
302    !--------- autres criteres
303
304    ! points poses
305    gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz)   !  Pas de calcul pour les points
306    gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz)   !  au fort beta
307
308
309
310    ! critere (lache) sur la temperature
311    do j = 2, ny-1
312       do i =2, nx-1 
313
314!   test s'il y a un point tempere dans les environs
315          test_Tx = any( ibase( i-1:i   , j-1:j+1 ).gt.1)
316          test_Ty = any( ibase( i-1:i+1 , j-1:j   ).gt.1)
317
318!   test s'il y a un point low drag dans les environs
319          test_beta_x = any( drag_centre( i-1:i   , j-1:j+1 ) .lt. beta_forc_fleuve )
320          test_beta_y = any( drag_centre( i-1:i+1 , j-1:j   ) .lt. beta_forc_fleuve )
321
322!   critere pour gzmx
323
324! Attention : les deux lignes suivants sont en commentaire pour voir l'effet
325
326!          gzmx(i,j) =  gzmx(i,j) .and. (test_Tx.or.test_beta_x)
327!          gzmy(i,j) =  gzmy(i,j) .and. (test_Ty.or.test_beta_y)
328
329       end do
330    end do
331
332
333    if (itracebug.eq.1)  call tracebug(' apres autres criteres ')
334    if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330) 
335
336    flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
337    flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
338
339    where ((.not.flotmx(:,:)).and.(.not.gzmx(:,:)))
340       betamx(:,:) = beta_limgz
341    end where
342
343    where ((.not.flotmy(:,:)).and.(.not.gzmy(:,:)))
344       betamy(:,:) = beta_limgz
345    end where
346
347        fleuvemx(:,:)=gzmx(:,:)
348        fleuvemy(:,:)=gzmy(:,:)
349   
350  end subroutine gzm_beta_prescr
351
352!______________________________________________________________________________________
353!>SUBROUTINE:   beta_min_value
354!! valeur mini que peut avoir beta (en Pa an m-1)
355!<
356
357  subroutine beta_min_value(val_betamin)
358
359  real :: val_betamin     ! valeur mini que peut avoir beta (en Pa an m-1)
360
361    drag_centre(:,:) = max(drag_centre(:,:),val_betamin)
362    drag_mx(:,:) = max(drag_mx(:,:),val_betamin)
363    drag_my(:,:) = max(drag_my(:,:),val_betamin)
364    where(flot(:,:))
365       drag_centre(:,:) = 0.
366    end where
367
368  end subroutine beta_min_value
369
370!>SUBROUTINE:   beta_c2stag 
371!! redistribute central values on staggered grid
372!<
373
374  subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre)              ! redistribute central values on staggered grid
375
376    implicit none
377    integer                                ::  nxx,nyy                ! dimension des tableaux
378    real, dimension(nxx,nyy), intent(out)  ::  R_mx                   ! valeur du beta sur maille mx
379    real, dimension(nxx,nyy), intent(out)  ::  R_my                   ! valeur du beta sur maille my
380    real, dimension(nxx,nyy), intent(in)   ::  R_centre               ! valeur du beta sur maille centree
381
382    do j=2,nyy 
383       do i=2,nxx
384          R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5
385          R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5
386       end do
387    end do
388
389  end subroutine beta_c2stag
390
391
392!-----------------------------------------------------------------------------------------
393!< subroutine : beta_iter_vitbil
394!< fait des iterations pour affiner la valeur de beta_centre pour s'ajuster sur vitbil
395!-----------------------------------------------------------------------------------------
396subroutine beta_iter_vitbil
397
398
399! calcule la vitesse centree venant du calcul direct
400
401  Umag_direct(:,:)    = (((uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))**2+ &        ! moyenne
402                          (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5
403
404  Uslmag_direct(:,:)  = (((ux(:,:,nz)+eoshift(ux(:,:,nz),shift=1,boundary=0.,dim=1))**2+ &        ! glissement
405                          (uy(:,:,nz)+eoshift(uy(:,:,nz),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5
406
407  Udefmag_direct(:,:) = max( (Umag_direct(:,:)- Uslmag_direct(:,:)) , 0.)                          ! deformation
408
409
410! guess du beta dans les zones où GRISLI a une vitesse de glissement nulle
411
412  call slope_surf
413  driving_stress(:,:) = rog * H(:,:) * slope(:,:)   !  driving stress
414
415
416  Uslid_vitbil(:,:)   = Umag_vitbil(:,:) - Udefmag_direct(:,:)  ! calcule la vitesse de glissement "bilan"
417
418
419! Pas de glissement dans GRISLI
420
421where ((Uslmag_direct(:,:).le.0.).and.(Uslid_vitbil(:,:).gt.1.e-3))     ! pas de glissement dans GRISLI
422  beta_centre(:,:) = driving_stress(:,:) / Uslid_vitbil(:,:)            ! mais du glissement "bilan"
423end where
424
425
426where ((Uslmag_direct(:,:).le.0.).and.(Uslid_vitbil(:,:).le.1.e-3))     ! ni l'un ni l'autre
427   Uslid_vitbil(:,:) = 0.
428   beta_centre(:,:)  = beta_limgz
429end where
430
431
432! Glissement dans GRISLI
433
434  where (Uslid_vitbil(:,:).le.1.e-3)    ! la vitesse de bilan est deja trop rapide
435     Uslid_vitbil(:,:) = 0.
436     beta_centre(:,:)  = beta_limgz
437
438  elsewhere (Uslmag_direct(:,:).gt.0.)
439     beta_centre(:,:) = beta_centre(:,:) * coef_iter_vitbil * Uslmag_direct(:,:) / Uslid_vitbil(:,:)
440
441  end where
442
443  where(flot(:,:))
444     beta_centre(:,:) = 0.
445  end where
446
447
448
449  beta_centre(:,:) = min(beta_centre(:,:),beta_limgz)
450 
451    call  beta_c2stag(nx,ny,betamx,betamy,beta_centre)  ! redistribue sur les mailles alternees
452    call  beta_min_value(beta_min)                      !  valeur mini que peut avoir beta (en Pa an m-1)
453 
454
455 debug_3D(:,:,59) = betamx(:,:)
456 debug_3D(:,:,60) = beta_centre(:,:)
457
458! estimateurs
459J_Umag  = 0.
460J_Udef   = 0.
461nb_umag = 0
462nb_udef = 0
463
464do j=2,ny-1
465   do i=2,nx-1
466
467      if (.not.flot(i,j)) then                       ! calcul sur les points poses
468         nb_umag = nb_umag + 1
469         J_Umag = (Umag_direct(i,j) - Umag_vitbil(i,j))**2. + J_Umag
470
471         if (Uslid_vitbil(i,j).le.1.e-3) then       ! calcul sur la partie deformation suelement
472            nb_udef = nb_udef + 1
473            J_Udef = (Umag_direct(i,j) - Umag_vitbil(i,j))**2. + J_Udef
474         end if
475      end if
476
477   end do
478end do
479      J_Umag = (J_Umag**0.5) / nb_umag
480      J_Udef = (J_Udef**0.5) / nb_udef
481
482
483drag_centre(:,:) = beta_centre(:,:)
484drag_mx(:,:)     = betamx(:,:)
485drag_my(:,:)     = betamy(:,:)
486
487i=319
488j=113
489write(266,'(2(i0,1x),3x,9(es12.4,1x))') nb_umag,nb_udef,time,J_Umag,J_Udef,Umag_direct(i,j),Umag_vitbil(i,j),Uslmag_direct(i,j),Uslid_vitbil(i,j),beta_centre(i,j),betamx(i,j),drag_centre(i,j)
490
491end subroutine beta_iter_vitbil
492!-----------------------------------------------------------------------------------------
493end module beta_iter_vitbil_mod
Note: See TracBrowser for help on using the repository browser.