source: branches/iLoveclim/SOURCES/Draggings_modules/dragging_beta_iter_vitbil_mod.f90 @ 187

Last change on this file since 187 was 187, checked in by aquiquet, 6 years ago

Grisli-iloveclim branch merged to trunk at revision 185
File size: 10.6 KB
Line 
1!< \file beta_iter_vitbil_mod.f90
2!< Module dragging 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_iter_vitbil_mod
5
6
7module dragging_beta_iter_vitbil_mod
8
9  use module3d_phy
10  use interface_input
11  use netcdf
12  use io_netcdf_grisli
13  use beta_iter_vitbil_mod
14
15  implicit none
16
17
18!~   real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite
19!~   real, dimension(nx,ny) :: Vcol_y           !<
20!~   real, dimension(nx,ny) :: Vsl_x            !<
21!~   real, dimension(nx,ny) :: Vsl_y            !<
22
23
24!~   real, dimension(nx,ny) :: drag_centre      !< beta on major node (average)
25!~   real :: beta_limgz                         !< when beta gt beta_limgz -> not gzmx
26!~   real :: beta_min                           !< minimum value of beta
27!~   real :: beta_mult                          !< coefficient of beta field
28!~   integer  :: ill,jll,nmoy
29
30!~   real :: maxi                               !< calcul correction dS a la grounding line
31!~   real :: mini
32
33!~   logical :: corr_def = .False.               !< for deformation correction (compatibility)
34
35
36!~   real, dimension(nx,ny) :: Umag_vitbil             !< vitesse de bilan
37!~   real, dimension(nx,ny) :: Uslid_vitbil            !< vitesse de bilan partie glissement
38!~   real, dimension(nx,ny) :: Umag_direct             !< vitesse moyenne centree calculee par GRISLI
39!~   real, dimension(nx,ny) :: Uslmag_direct           !< vitesse glissement centree calculee par GRISLI
40!~   real, dimension(nx,ny) :: Udefmag_direct          !< vitesse deformation centree calculee par GRISLI
41!~   real, dimension(nx,ny) :: driving_stress          !< driving stress
42
43!~   character(len=100)     :: Umag_bil_file           !< fichier qui contient les donnees
44!~   real                   :: time_iter               !< temps de debut des iterations
45!~   real                   :: time_iter_end           !< temps de fin des iterations
46!~   real                   :: time_reiter             !< nbr annees entre 2 iterations calcul beta
47!~   integer                :: nb_iter_vitbil          !< nombre d'iteratiosn avant de changer de pas de temps
48!~   real                   :: coef_iter_vitbil        !< coefficient <= 1 (rapport vitesses)
49!~   real                   :: J_Umag                  !< estimateur ecart entre vitesses
50!~   real                   :: J_Udef                  !< estimateur vitesse regions sans glissement
51!~   integer                :: nb_umag                 !< nombre de points calcul J_Umag
52!~   integer                :: nb_Udef                 !< nombre de points calcul J_Umag
53
54
55contains
56
57  !-----------------------------------------------------------------------------------------
58  !> SUBROUTINE: init_dragging
59  !! Cette routine fait l'initialisation du dragging.
60  !<
61  subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
62    ! nouvelle version : lit les fichiers nc
63    implicit none
64    character(len=100) :: beta_c_file               !  beta on centered grid
65    character(len=100) :: betamx_file               !  beta on staggered grid mx
66    character(len=100) :: betamy_file               !  beta on staggered grid mx
67    real*8, dimension(:,:),   pointer  :: tab => null() !< tableau 2d real ecrit dans le fichier pour lect ncdf
68
69    namelist/beta_prescr/beta_c_file,beta_limgz,beta_min,beta_mult
70
71    if (itracebug.eq.1)  call tracebug('beta_iter_vitbil_mod subroutine init_dragging')
72
73    inv_beta=1
74
75    ! lecture des parametres du run                   
76    !--------------------------------------------------------------------
77
78    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
79428 format(A)
80    read(num_param,beta_prescr)
81
82    write(num_rep_42,428) '!___________________________________________________________' 
83    write(num_rep_42,428) '!  read beta on centered grid'
84    write(num_rep_42,beta_prescr)                   
85    write(num_rep_42,428) '!  beta_file : nom des fichier qui contiennent les beta_c'
86    write(num_rep_42,428) '!  above beta_limgz, gzmx is false'
87    write(num_rep_42,428) '!  if grounded, beta_min minimum value '
88    write(num_rep_42,428) '!  beta_mult : coefficient just after reading' 
89    write(num_rep_42,*)
90    write(num_rep_42,428) '!___________________________________________________________' 
91
92
93    beta_c_file = trim(dirnameinp)//trim(beta_c_file)
94    betamx_file = trim(dirnameinp)//trim(betamx_file)
95    betamy_file = trim(dirnameinp)//trim(betamy_file)
96    betamax = beta_limgz
97    betamax_2d(:,:) = betamax
98
99!    call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'')
100    write(6,*)  beta_c_file
101
102    call Read_Ncdf_var('z',beta_c_file,tab)
103    drag_centre(:,:) = tab(:,:)
104
105
106    ! multiplication par le coefficient beta_mult
107
108    drag_centre(:,:) = drag_centre(:,:)*beta_mult
109
110    where (mk_init(:,:).eq.-2)                                             ! iles a probleme
111       drag_centre(:,:) = 2.* abs(beta_limgz)
112    end where
113
114    where(flot(:,:))
115       drag_centre(:,:) = 0.
116    end where
117
118
119    call  beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)         ! redistribute on staggered grid
120
121    if  (beta_limgz.gt.0.) then
122       beta_centre(:,:) =  drag_centre(:,:)
123       betamx(:,:)      =  drag_mx(:,:)
124       betamy(:,:)      =  drag_my(:,:)
125
126    else if (beta_limgz.lt.0.) then                              ! attention sans doute pour plastic
127
128       beta_centre(:,:) = - beta_limgz
129       betamx(:,:)      = - beta_limgz
130       betamy(:,:)      = - beta_limgz
131       drag_centre(:,:) = - beta_limgz 
132       drag_mx(:,:)     = - beta_limgz
133       drag_my(:,:)     = - beta_limgz 
134       beta_limgz = 1.
135    end if
136
137
138
139    call beta_min_value(beta_min)      !  valeur mini que peut avoir beta (en Pa an m-1)
140    ! on peut envisager une valeur ~ 10
141    ! rappel : beta doit être positif.
142
143
144
145    !    call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)                  ! centered distributed
146    !                                                                         ! on staggered
147    !
148
149
150    beta_centre(:,:) =  drag_centre(:,:)                                 ! surtout pour sorties
151
152
153
154    where (drag_mx(:,:).ge.beta_limgz)
155       mstream_mx(:,:) = 0
156       betamx(:,:)     = beta_limgz
157       drag_mx(:,:)    = beta_limgz
158    elsewhere
159       mstream_mx(:,:) = 1
160       betamx(:,:)     = drag_mx(:,:)
161    endwhere
162
163    where (drag_my(:,:).ge.beta_limgz)
164       mstream_my(:,:) = 0
165       betamy(:,:)     = beta_limgz
166       drag_my(:,:)    = beta_limgz
167    elsewhere
168       mstream_my(:,:) = 1
169       betamy(:,:)     = drag_my(:,:)
170    endwhere
171
172    if (itracebug.eq.1)  call tracebug(' Prescr_beta  apres lecture')
173    mstream_mx(:,:) = 0
174    mstream_my(:,:) = 0
175
176
177    call gzm_beta_prescr
178
179    call init_beta_iter_vitbil
180    return
181
182  end subroutine init_dragging
183
184  !________________________________________________________________________________
185
186  !> SUBROUTINE: dragging
187  !! Defini la localisation des streams et le frottement basal
188  !<
189  !-------------------------------------------------------------------------
190  subroutine dragging   ! defini la localisation des streams et le frottement basal
191
192    implicit none
193
194    if (abs(time-time_iter_end).lt.dtmin) then
195      time_iter = time_iter + time_reiter
196      time_iter_end = time_iter_end + time_reiter
197      print*,'debug dragging beta_iter_vitbil_mod time_iter=', time_iter, time_iter_end, time
198    endif
199
200    if (itracebug.eq.1)  call tracebug(' beta_prescr         subroutine dragging')
201
202
203    betamx(:,:)=drag_mx(:,:)
204    betamy(:,:)=drag_my(:,:)
205
206
207    call  gzm_beta_prescr        ! determine gz et flgz et met a 0 le beta des points flottants
208
209
210    return
211  end subroutine dragging
212
213  !----------------------------------------------------------------------------------
214  !>SUBROUTINE: gzm_beta_prescr
215  !! Calcul de gzmx
216  !<
217
218  subroutine gzm_beta_prescr     ! attribue flgzmx et gzmx (et my)
219   
220    logical ::  test_Tx          ! test if there is a basal melting point in the surrounding
221    logical ::  test_Ty          ! test if there is a basal melting point in the surrounding
222    logical ::  test_beta_x      ! test if there is a low drag point in the surrounding
223    logical ::  test_beta_y      ! test if there is a low drag point in the surrounding
224
225    real    ::  beta_forc_fleuve ! below this value -> ice stream
226
227!    beta_forc_fleuve = 5.e3
228    beta_forc_fleuve = beta_limgz
229
230    ! calcul de gzmx
231
232
233    gzmx(:,:)=.true.
234    gzmy(:,:)=.true.
235    flgzmx(:,:)=.false.
236    flgzmy(:,:)=.false.
237
238
239    ! points flottants : flgzmx mais pas gzmx
240    do j=2,ny
241       do i=2,nx
242          if (flot(i,j).and.(flot(i-1,j))) then
243             flotmx(i,j)=.true.
244             flgzmx(i,j)=.true.
245             gzmx(i,j)=.false.
246             betamx(i,j)=0.
247          end if
248          if (flot(i,j).and.(flot(i,j-1))) then
249             flotmy(i,j)=.true.
250             flgzmy(i,j)=.true.
251             gzmy(i,j)=.false.
252             betamy(i,j)=0.
253          end if
254       end do
255    end do
256
257    if (itracebug.eq.1)  call tracebug(' apres criteres flot')
258
259    if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,127) 
260
261    !--------- autres criteres
262
263    ! points poses
264    gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz)   !  Pas de calcul pour les points
265    gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz)   !  au fort beta
266
267
268
269    ! critere (lache) sur la temperature
270    do j = 2, ny-1
271       do i =2, nx-1 
272
273!   test s'il y a un point tempere dans les environs
274          test_Tx = any( ibase( i-1:i   , j-1:j+1 ).gt.1)
275          test_Ty = any( ibase( i-1:i+1 , j-1:j   ).gt.1)
276
277!   test s'il y a un point low drag dans les environs
278          test_beta_x = any( drag_centre( i-1:i   , j-1:j+1 ) .lt. beta_forc_fleuve )
279          test_beta_y = any( drag_centre( i-1:i+1 , j-1:j   ) .lt. beta_forc_fleuve )
280
281!   critere pour gzmx
282
283! Attention : les deux lignes suivants sont en commentaire pour voir l'effet
284
285!          gzmx(i,j) =  gzmx(i,j) .and. (test_Tx.or.test_beta_x)
286!          gzmy(i,j) =  gzmy(i,j) .and. (test_Ty.or.test_beta_y)
287
288       end do
289    end do
290
291
292    if (itracebug.eq.1)  call tracebug(' apres autres criteres ')
293    if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,127) 
294
295    flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
296    flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
297
298    where ((.not.flotmx(:,:)).and.(.not.gzmx(:,:)))
299       betamx(:,:) = beta_limgz
300    end where
301
302    where ((.not.flotmy(:,:)).and.(.not.gzmy(:,:)))
303       betamy(:,:) = beta_limgz
304    end where
305
306        fleuvemx(:,:)=gzmx(:,:)
307        fleuvemy(:,:)=gzmy(:,:)
308   
309  end subroutine gzm_beta_prescr
310
311
312!-----------------------------------------------------------------------------------------
313end module dragging_beta_iter_vitbil_mod
Note: See TracBrowser for help on using the repository browser.