source: trunk/SOURCES/dragging_neff_slope_mod.f90 @ 66

Last change on this file since 66 was 66, checked in by aquiquet, 8 years ago

Dragging neff parameters in param_list && distinction between actual ice streams and SSA velocity used as Ub in SIA (slow SSA zones) in netcdf outputs

File size: 14.1 KB
Line 
1!> \file dragging_neff_contmaj_mod.f90
2!! Defini les zones de stream avec 3 criteres
3!<
4
5!> \namespace  dragging_neff_contmaj
6!! Defini les zones de stream avec 3 criteres
7!! \author ...
8!! \date ...
9!! @note Les trois criteres sont:
10!! @note   * un critere sur la pression effective
11!! @note   * un critere de continuite depuis la cote
12!! @note   * un critere sur la courbure du socle (si negatif, vallees)
13!! @note Used module
14!! @note   - use module3D_phy
15!! @note   - use param_phy_mod
16!<
17
18module dragging_neff_slope
19
20! Defini les zones de stream avec :
21! * un critere sur la pression effective
22! * un critere de continuite depuis la cote
23! * un critere sur la courbure du socle (si negatif, vallees)
24
25use module3d_phy
26use param_phy_mod
27use interface_input
28use io_netcdf_grisli
29
30implicit none
31
32logical,dimension(nx,ny) :: fleuve
33logical,dimension(nx,ny) :: cote
34logical,dimension(nx,ny) :: slowssamx ! not actual stream, but ssa used as Ub
35logical,dimension(nx,ny) :: slowssamy ! not actual stream, but ssa used as Ub
36logical,dimension(nx,ny) :: slowssa   ! not actual stream, but ssa used as Ub
37
38real,dimension(nx,ny) :: hires_slope ! slope comupted on a high resolution topography
39character(len=100) :: slope_fich     ! fichier grille
40character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat
41
42real :: valmax
43integer :: imax,jmax
44integer :: i_moins1,i_plus1,j_moins1,j_plus1
45integer :: lmax=20
46integer :: idep,jdep,iloc,jloc
47
48
49real :: tostick   ! pour la glace posee
50real :: betamin   ! betamin : (Pa) frottement mini sous les streams
51real :: tob_ile    ! pour les iles
52real :: cry_lim=50.  ! courbure limite pour le suivi des fleuves
53
54real,dimension(nx,ny) :: neff ! pression effective noeuds majeurs
55real :: seuil_neff    ! seuil sur la pression effective pour avoir glissement en zone sediment
56real :: coef_gz       ! coef frottement zones stream std  (10)
57real :: coef_ile      ! coef frottement zones iles   (0.1)
58
59real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
60real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
61real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
62real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
63logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta
64
65contains
66!-------------------------------------------------------------------------------
67
68!> SUBROUTINE: init_dragging
69!! Cette routine fait l'initialisation du dragging.
70!>
71subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
72
73implicit none
74
75! local variables, defining the rugosity-enhanced flow
76real :: expo_slope
77real :: pente_min, pente_max
78
79namelist/drag_neff_slope/cf,betamax,betamin,toblim,tostick,seuil_neff,coef_gz,coef_ile,slope_fich,expo_slope,pente_min,pente_max
80
81if (itracebug.eq.1)  call tracebug(' dragging avec neff et slope')
82
83
84! formats pour les ecritures dans 42
85428 format(A)
86
87! lecture des parametres du run                      block drag neff
88!--------------------------------------------------------------------
89
90rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
91read(num_param,drag_neff_slope)
92
93write(num_rep_42,428)'!___________________________________________________________' 
94write(num_rep_42,428) '&drag_neff_slope          ! nom du bloc dragging neff slope'
95write(num_rep_42,*)
96write(num_rep_42,*) 'cf               = ',cf
97write(num_rep_42,*) 'betamax          = ', betamax
98write(num_rep_42,*) 'betamin          = ', betamin
99write(num_rep_42,*) 'toblim           = ', toblim
100write(num_rep_42,*) 'tostick          = ', tostick
101write(num_rep_42,*) 'seuil_neff       = ', seuil_neff
102write(num_rep_42,*) 'coef_gz          = ', coef_gz
103write(num_rep_42,*) 'coef_ile         = ', coef_ile
104write(num_rep_42,'(A,A)') 'slope_fich = ', slope_fich
105write(num_rep_42,*) 'expo_slope       = ', expo_slope
106write(num_rep_42,*) 'pente_min        = ', pente_min
107write(num_rep_42,*) 'pente_max        = ', pente_max
108write(num_rep_42,*)'/'                           
109write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
110write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams '
111write(num_rep_42,428) '! betamin : (Pa) frottement mini sous les streams '
112write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
113write(num_rep_42,428) '! tostick : (Pa)  pour les points non flgzmx '
114write(num_rep_42,428) '! seuil_neff (Pa) seuil sur la pression effective pour avoir glissement'
115write(num_rep_42,428) '! coef_gz : coef frottement zones stream std'
116write(num_rep_42,428) '! coef_ile : coef frottement zones iles'
117write(num_rep_42,428) '! slope_fich : fichier pente bedrock'
118write(num_rep_42,428) '! expo_slope : exposant fonction pente'
119write(num_rep_42,428) '! pente_min : no flow enhancement below this'
120write(num_rep_42,428) '! pente_max : maximum flow enhancement above this'
121write(num_rep_42,*)
122
123!tostick=1.e5   ! valeurs pour les points non flgzmx
124tob_ile=betamax/2.
125moteurmax=toblim
126
127! betamax_2d depends on sub-grid slopes.
128slope_fich=trim(dirnameinp)//trim(slope_fich)
129call lect_input(1,'slope',1,hires_slope(:,:),slope_fich,file_ncdf)
130
131! from slopes, we create an index between 0 & 1
132! 1 is very mountainous, 0 is flat
133hires_slope(:,:)=(max(min(hires_slope(:,:),pente_max),pente_min)-pente_min)/(pente_max-pente_min)
134! now we compute the actual betamax used by the remplimat routines
135! /!\ the slope is used to modify the beta where we have a temperate base,
136! but NO ice stream... -> Slow SSA zone (SSA used to compute Ub)
137betamax_2d(:,:)=max ( tostick * (1. - (1 - betamax / tostick ) * hires_slope(:,:)**(1./expo_slope)) , betamax )
138!do j=1,ny
139!   do i=1,nx
140!      write(18745,*) betamax_2d (i,j)
141!   enddo
142!enddo
143
144!-------------------------------------------------------------------
145! masque stream
146      mstream_mx(:,:)=1
147      mstream_my(:,:)=1
148
149
150! coefficient permettant de modifier le basal drag.
151      drag_mx(:,:)=1.
152      drag_my(:,:)=1.
153
154where (socle_cry(:,:).lt.cry_lim)
155   debug_3D(:,:,25)=1
156elsewhere
157   debug_3D(:,:,25)=0
158endwhere
159
160   
161      return
162      end subroutine init_dragging
163!________________________________________________________________________________
164
165!> SUBROUTINE: dragging
166!! Defini la localisation des streams et le frottement basal
167!>
168!-------------------------------------------------------------------------
169subroutine dragging   ! defini la localisation des streams et le frottement basal
170
171
172!         les iles n'ont pas de condition neff mais ont la condition toblim
173!         (ce bloc etait avant dans flottab)
174
175
176do j=2,ny
177   do i=2,nx
178      ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim)
179      ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim)
180   end do
181end do
182
183! le calcul des fleuves se fait sur les mailles majeures en partant de la cote
184
185
186
187
188! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
189
190! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max
191! tant que ce point est superieur a hwatermax.
192
193! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
194!             et les fleuves n'ont une épaisseur que de 1 noeud.
195
196fleuvemx(:,:)=.false.
197fleuvemy(:,:)=.false.
198fleuve(:,:)=.false.
199cote(:,:)=.false.
200
201! detection des cotes
202do  j=2,ny-1 
203   do i=2,nx-1 
204      if ((.not.flot(i,j)).and. & 
205      ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then
206         cote(i,j)=.true.
207      endif
208   end do
209end do
210
211! calcul de neff (pression effective sur noeuds majeurs)
212do j=1,ny-1
213   do i=1,nx-1
214        neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4
215   enddo
216enddo
217!aurel, for the borders:
218neff(:,ny)=1.e5
219neff(nx,:)=1.e5
220
221!!$
222!!$fleuve_maj: do j=2,ny-1 
223!!$ifleuve:   do i=2,nx-1                     
224!!$
225!!$cote_detect :  if (cote(i,j)) then
226!!$         idep=i
227!!$         jdep=j
228!!$
229!!$         if (socle_cry(i,j).lt.0.) then    ! dans une vallee
230!!$            fleuve(i,j)=.true.
231!!$         else
232!!$            cote(i,j)=.false.
233!!$            cycle ifleuve
234!!$         endif
235!!$
236!!$suit : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
237!!$           i_moins1=max(idep-1,2)
238!!$           j_moins1=max(jdep-1,1)
239!!$           i_plus1=min(idep+1,nx)
240!!$           j_plus1=min(jdep+1,ny)
241!!$           
242!!$! recherche du max en suivant le socle le plus profond
243!!$! * en excluant les points flottants
244!!$! * et ceux qui sont deja tagges fleuves
245!!$
246!!$           valmax=1000.
247!!$
248!!$           do jloc=j_moins1,j_plus1
249!!$              do iloc=i_moins1,i_plus1
250!!$
251!!$                 if ((B(iloc,jloc).lt.valmax)      &
252!!$                      .and.(.not.flot(iloc,jloc))  &
253!!$                      .and.(.not.fleuve(iloc,jloc)).and.(socle_cry(iloc,jloc).lt.cry_lim)) then
254!!$                    imax=iloc
255!!$                    jmax=jloc
256!!$                    valmax=B(iloc,jloc)
257!!$                 endif
258!!$              end do
259!!$           end do
260!!$
261!!$         if ((hwater(imax,jmax).gt.hwatstream).and.(socle_cry(i,j).lt.cry_lim)) then
262!!$            fleuve(imax,jmax)=.true.
263!!$            idep=imax
264!!$            jdep=jmax
265!!$         else
266!!$            fleuve(imax,jmax)=.false.
267!!$            exit suit
268!!$         end if
269!!$
270!!$      end do suit
271!!$
272!!$   end if cote_detect
273!!$
274!!$end do ifleuve
275!!$end do fleuve_maj
276
277! aurel, we add the neff threshold:
278where ((neff(:,:).le.seuil_neff).and.(.not.flot(:,:)).and.(H(:,:).gt.1.)) fleuve(:,:)=.true.
279
280do j=1,ny-1
281   do i=1,nx-1
282      if (fleuve(i,j)) then
283         fleuvemx(i,j)=.true.
284         fleuvemx(i+1,j)=.true.
285         fleuvemy(i,j)=.true.
286         fleuvemy(i,j+1)=.true.
287      end if
288   end do
289end do
290
291! we look for the warm based points that will not be treated as stream (ub from SSA):
292slowssa(:,:)=.false.
293slowssamx(:,:)=.false.
294slowssamy(:,:)=.false.
295do  j=1,ny
296   do i=1,nx
297      !if ((not(flot(i,j))).and.(hwater(i,j).gt.0.1)) slowssa(i,j)=.true.
298      if ((not(flot(i,j))).and.(ibase(i,j).ne.1).and.(H(i,j).gt.1.)) slowssa(i,j)=.true.
299   end do
300end do
301do j=1,ny-1
302   do i=1,nx-1
303      if (slowssa(i,j)) then
304         slowssamx(i,j)=.true.
305         slowssamx(i+1,j)=.true.
306         slowssamy(i,j)=.true.
307         slowssamy(i,j+1)=.true.
308      end if
309   end do
310end do
311
312do j=1,ny
313   do i=1,nx
314
315      ! the actual streams and the warm based points are gzmx now:
316      if ( ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) .or. ((.not.ilemx(i,j)).and.(slowssamx(i,j))) ) gzmx(i,j)=.true. 
317
318     
319! calcul du frottement basal (ce bloc etait avant dans neffect)
320
321      if (cotemx(i,j)) then                        ! point cotier
322         betamx(i,j)=cf*neffmx(i,j) 
323         betamx(i,j)=min(betamx(i,j),betamax)
324
325      else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then ! tous les points SSA
326
327         if (fleuvemx(i,j)) then                 ! the actual streams
328            betamx(i,j)=cf*neffmx(i,j)*coef_gz
329            betamx(i,j)=min(betamx(i,j),betamax)
330            betamx(i,j)=max(betamx(i,j),betamin)
331
332            if (cf*neffmx(i,j).gt.betamax*2.) then  ! a stream that's becoming grounded...
333               if (slowssamx(i,j)) then
334                  betamx(i,j)=betamax_2d(i,j)
335               else
336                  gzmx(i,j)=.false.
337                  betamx(i,j)=tostick
338               endif
339            endif
340
341         else                    ! not an actual stream, SSA is used to compute Ub
342            betamx(i,j)=betamax_2d(i,j)
343         endif
344         
345      else if (ilemx(i,j)) then
346         betamx(i,j)=cf*neffmx(i,j)*coef_ile 
347         betamx(i,j)=min(betamx(i,j),tob_ile)       
348      else if (flotmx(i,j)) then     ! flottant ou ile
349         betamx(i,j)=0.
350      else                                         ! grounded, SIA
351         betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
352      endif
353
354   end do
355end do
356
357do j=1,ny
358   do i=1,nx
359
360      ! the actual streams and the warm based points are gzmx now:
361      if ( ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) .or. ((.not.ilemy(i,j)).and.(slowssamy(i,j))) ) gzmy(i,j)=.true. 
362
363     
364! calcul du frottement basal (ce bloc etait avant dans neffect)
365
366      if (cotemy(i,j)) then                        ! point cotier
367         betamy(i,j)=cf*neffmy(i,j) 
368         betamy(i,j)=min(betamy(i,j),betamax)
369
370      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then ! tous les points SSA
371
372         if (fleuvemy(i,j)) then                 ! the actual streams
373            betamy(i,j)=cf*neffmy(i,j)*coef_gz
374            betamy(i,j)=min(betamy(i,j),betamax)
375            betamy(i,j)=max(betamy(i,j),betamin)
376
377            if (cf*neffmy(i,j).gt.betamax*2.) then  ! a stream that's becoming grounded...
378               if (slowssamy(i,j)) then
379                  betamy(i,j)=betamax_2d(i,j)
380               else
381                  gzmy(i,j)=.false.
382                  betamy(i,j)=tostick
383               endif
384            endif
385
386         else                    ! not an actual stream, SSA is used to compute Ub
387            betamy(i,j)=betamax_2d(i,j)
388         endif
389         
390      else if (ilemy(i,j)) then
391         betamy(i,j)=cf*neffmy(i,j)*coef_ile 
392         betamy(i,j)=min(betamy(i,j),tob_ile)       
393      else if (flotmy(i,j)) then     ! flottant ou ile
394         betamy(i,j)=0.
395      else                                         ! grounded, SIA
396         betamy(i,j)=tostick                       ! frottement glace posee (1 bar)
397      endif
398
399   end do
400end do
401
402
403do j=2,ny-1
404   do i=2,nx-1
405      beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) &
406          + (betamy(i,j)+betamy(i,j+1)))*0.25
407   end do
408end do
409
410
411where (fleuve(:,:))
412   debug_3D(:,:,1)=1
413elsewhere (flot(:,:))
414   debug_3D(:,:,1)=2
415elsewhere
416   debug_3D(:,:,1)=0
417endwhere
418
419where (cote(:,:))
420   debug_3D(:,:,2)=1
421elsewhere
422   debug_3D(:,:,2)=0
423endwhere
424
425where (fleuvemx(:,:))
426   debug_3D(:,:,3)=1
427elsewhere
428   debug_3D(:,:,3)=0
429endwhere
430
431where (flotmx(:,:))
432   debug_3D(:,:,3)=-1
433endwhere
434
435!_____________________
436
437
438where (fleuvemy(:,:))
439   debug_3D(:,:,4)=1
440elsewhere
441   debug_3D(:,:,4)=0
442endwhere
443
444where (flotmy(:,:))
445   debug_3D(:,:,4)=-1
446end where
447
448debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) 
449debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) 
450
451return
452end subroutine dragging
453
454end module dragging_neff_slope
455
Note: See TracBrowser for help on using the repository browser.