source: trunk/SOURCES/furst_schoof_mod.f90 @ 271

Last change on this file since 271 was 271, checked in by aquiquet, 5 years ago

Further computations of mass conservation in double precision for small timesteps
File size: 26.7 KB
Line 
1! Module Furst Schoof
2
3module furst_schoof_mod
4
5use module3d_phy
6use deformation_mod_2lois
7
8implicit none
9
10real :: frot_coef ! afq: fixed in Tsai but = beta in Schoof (for m_weert = 1)
11real :: alpha_flot
12integer :: gr_select
13real,dimension(nx,ny) :: back_force_x
14real,dimension(nx,ny) :: back_force_y
15real, parameter :: m_weert = 1.0 ! afq: is this defined elsewhere in the model???
16real, parameter :: inv_mweert = 1./m_weert
17
18contains
19
20subroutine init_furst_schoof
21
22namelist/furst_schoof/frot_coef,gr_select
23
24rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
25read(num_param,furst_schoof)
26write(num_rep_42,'(A)')'!___________________________________________________________' 
27write(num_rep_42,'(A)') '&furst_schoof                                ! nom du bloc '
28write(num_rep_42,*)
29write(num_rep_42,*)'frot_coef  = ',frot_coef
30write(num_rep_42,*)'gr_select  = ',gr_select
31write(num_rep_42,*)'/'
32write(num_rep_42,*)'! frot_coef : solid friction law coef 0.6 in GMD 2018'
33write(num_rep_42,*)'! gr_select = 1 : Tsai , 2 : Schoof'
34
35!  frot = 0.035 ! f in the solid friction law, 0.035 in Ritz et al. 2001       
36
37  ! choix Tsai (loi de Coulomb) ou Schoof (loi de Weertman)
38  alpha_flot= ro/row
39  back_force_x(:,:)=0.1
40  back_force_y(:,:)=0.1
41 
42end subroutine init_furst_schoof 
43 
44subroutine interpol_glflux
45
46  implicit none
47  real :: xpos,ypos
48  real :: Hgl
49  !  integer :: igl
50
51  real :: phi_im2, phi_im1, phi_i, phi_ip1, phi_ip2, phi_ip3
52  real :: phi_jm2, phi_jm1, phi_j, phi_jp1, phi_jp2, phi_jp3
53  real,dimension(nx,ny) :: phi_xgl, phi_ygl ! flux gr line sur maille stag
54  integer,dimension(nx,ny) :: countx, county   ! how often do we modify ux/uy
55  real :: phi_prescr
56  real :: toutpetit = 1e-6
57  real :: denom, prodscal
58
59  real :: bfx, bfy
60
61  !debug
62  real,dimension(nx,ny) :: xpos_tab=9999.
63  real,dimension(nx,ny) :: ypos_tab=9999.
64  real,dimension(nx,ny) :: Hglx_tab=9999.
65  real,dimension(nx,ny) :: Hgly_tab=9999.
66  real,dimension(nx,ny) :: Hmxloc, Hmyloc ! pour eviter division par 0
67  real,dimension(nx,ny) :: phi_prescr_tabx, phi_prescr_taby
68
69  real,dimension(nx,ny) :: uxbartemp, uybartemp
70 
71  real,dimension(nx,ny) :: archimtab
72 
73  real,dimension(nx,ny) :: imx_diag_in
74
75  phi_prescr_tabx=0.
76  phi_prescr_taby=0.
77
78  Hmxloc(:,:)= max(Hmx(:,:),1.)
79  Hmyloc(:,:)=max(Hmy(:,:),1.)
80
81  uxbartemp(:,:)=0.
82  uybartemp(:,:)=0.
83  countx(:,:)=0.
84  county(:,:)=0.
85  gr_line_schoof(:,:)=0
86
87  archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel_2d(:,:)
88
89  imx_diag_in(:,:) = imx_diag(:,:)
90 
91  ! detection des noeuds grounding line
92  do j= 3,ny-3
93     do i=3,nx-3
94        !archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel
95        if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then    ! grounded and with ice
96           !           print*,'schoof gr line ij',i,j
97           ! selon x (indice i)
98           
99           if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) then
100           !if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then    ! grounding line between i-1,i    CAS WEST (ecoulement vers grille West)
101
102              call  calc_xgl4schoof(-dx,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i-1,j)),real(Bsoc(i-1,j)),sealevel_2d(i-1,j),xpos,Hgl)
103              xpos_tab(i,j)=xpos
104              Hglx_tab(i,j)=Hgl
105             
106              ! afq: the back force is on the staggered grid, the GL can be either West or East to this point.
107              if (xpos .lt. -dx/2.) then
108                 bfx = back_force_x(i,j)
109              else
110                 bfx = back_force_x(i+1,j)
111              endif
112              if (gr_select.eq.1) then ! flux de Tsai
113                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr)
114              else if (gr_select.eq.2) then ! flux de Schoof
115                 ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point.
116                 if (xpos .lt. -dx/2.) then
117                    frot_coef = betamx(i,j)
118                 else
119                    frot_coef = betamx(i+1,j)
120                 endif
121                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr)
122              else
123                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
124                 stop
125              endif
126              phi_prescr=-phi_prescr ! Doit être negatif dans le cas West
127              phi_prescr_tabx(i,j)=phi_prescr
128
129              if (xpos .lt. -dx/2.) then  !  GL a west du i staggered,                          o centre, x stag
130
131                 !  grille    ! .....x......o......x......o...G..x......O......x......o......x......o
132                 !            !            i-2           i-1            i            i+1           i+2
133                 !flux        !     in            out         G out            in
134
135
136!afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i-1,j))) then
137                   
138                 ! extrapolation pour avoir uxbar(i-1,j) qui sera ensuite prescrit dans call rempli_L2
139
140                 phi_im2 = Uxbar(i-2,j) * Hmxloc(i-2,j)
141                 phi_xgl(i-1,j)  = phi_im2 + dx * (phi_prescr-phi_im2)/(2.5 * dx + xpos)
142                 uxbartemp(i-1,j)   = uxbartemp(i-1,j)+ phi_xgl(i-1,j) / Hmxloc(i-1,j) ! attention division par 0 possible ?
143                 countx   (i-1,j)   = countx(i-1,j)+1
144                 imx_diag (i-1,j)   = 1
145                 gr_line_schoof(i-1,j) = 1
146
147                 ! extrapolation pour avoir uxbar(i,j)
148
149                 phi_ip1 = Uxbar(i+1,j) * Hmxloc(i+1,j)
150                 phi_xgl(i,j) = phi_ip1 + dx * (phi_prescr - phi_ip1)/(0.5 * dx - xpos)
151                 !                 print*,'uxbar(i,j)=',uxbar(i,j)
152                 uxbartemp(i,j)   = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j)
153                 countx   (i,j)   = countx(i,j)+1
154                 imx_diag (i,j)   = 1
155                 gr_line_schoof(i,j) = 1
156                 !                 print*,'schoof uxbar(i,j)=',uxbar(i,j)
157                 !                 print*,'Hgl',Hgl
158                 !                 print*,'phi_prescr',phi_prescr, phi_prescr/Hgl
159                 !                 print*
160
161!afq --                 end if
162             
163              else            !  GL a Est du i staggered,                          o centre, x stag
164
165                 !  grille    ! .....x......o......x......o......x...G..O......x......o......x......o
166                 !            i-2           i-1            i            i+1           i+2
167                 !flux        !     in             in           out  G         out           in
168
169!afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then
170
171                 ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
172
173                 phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j)
174                 phi_xgl(i,j)  = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos)
175                 uxbartemp(i,j)   = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ?
176                 countx   (i,j)   = countx(i,j)+1
177                 imx_diag (i,j) = 1
178                 gr_line_schoof(i,j) = 1
179
180                 ! extrapolation pour avoir uxbar(i+1,j)
181
182                 phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j)
183                 phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5 * dx  - xpos)
184                 uxbartemp(i+1,j)   = uxbartemp(i+1,j) +  phi_xgl(i+1,j) / Hmxloc(i+1,j)
185                 countx   (i+1,j)   = countx(i+1,j)+1
186                 imx_diag (i+1,j)   = 1
187                 gr_line_schoof(i+1,j) = 1
188
189!afq --                 end if
190
191              end if
192
193           else if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) then
194           !else if (flot(i+1,j).and.Uxbar(i+1,j).GT.0.) then    ! grounding line between i,i+1    CAS EST (ecoulement vers grille Est)
195
196              call  calc_xgl4schoof(dx,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i+1,j)),real(Bsoc(i+1,j)),sealevel_2d(i+1,j),xpos,Hgl)
197              xpos_tab(i,j)=xpos
198              Hglx_tab(i,j)=Hgl
199             
200              ! afq: the back force is on the staggered grid, the GL can be either West or East to this point.
201              if (xpos .lt. dx/2.) then
202                 bfx = back_force_x(i,j)
203              else
204                 bfx = back_force_x(i+1,j)
205              endif
206              if (gr_select.eq.1) then ! flux de Tsai
207                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr)
208              else if (gr_select.eq.2) then ! flux de Schoof
209                 ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point.
210                 if (xpos .lt. dx/2.) then
211                    frot_coef = betamx(i,j)
212                 else
213                    frot_coef = betamx(i+1,j)
214                 endif
215                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr)
216              else
217                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
218                 stop
219              endif
220              phi_prescr_tabx(i,j)=phi_prescr
221              if (xpos .lt. dx/2.) then  !  GL a west du i staggered,                          o centre, x stag
222
223                 !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o
224                 !            !            i-2           i-1            i            i+1           i+2
225                 !flux        !                    in           out        G  out            in
226
227!afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then
228                   
229                 ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
230
231                 phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j)
232                 phi_xgl(i,j)  = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos)
233                 uxbartemp(i,j)   =  uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ?
234                 countx   (i,j)   = countx(i,j)+1
235                 imx_diag (i,j) = 1
236                 gr_line_schoof(i,j) = 1
237
238                 ! extrapolation pour avoir uxbar(i+1,j)
239
240                 phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j)
241                 phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5* dx - xpos)
242                 uxbartemp(i+1,j)   = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j)
243                 countx   (i+1,j)   = countx(i+1,j)+1
244                 imx_diag (i+1,j) = 1
245                 gr_line_schoof(i+1,j) = 1
246
247!afq --                 end if
248             
249              else   !  GL a west du i staggered,                          o centre, x stag
250
251                 !  grille    ! ......x......o......x......O......x...G..o......x......o......x.....o
252                 !            !             i-1            i            i+1           i+2          i+3
253                 !flux        !                     in           out  G        out            in
254
255!afq --                 if (.not.(ilemx(i+1,j)).and..not.(ilemx(i+2,j))) then
256                   
257                 ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
258
259                 phi_i = Uxbar(i,j) * Hmxloc(i,j)
260                 phi_xgl(i+1,j)  = phi_i + dx * (phi_prescr-phi_i)/(0.5 * dx + xpos)
261                 uxbartemp(i+1,j)   = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j) ! attention division par 0 possible ?
262                 countx   (i+1,j)   = countx(i+1,j)+1
263                 imx_diag (i+1,j) = 1
264                 gr_line_schoof(i+1,j) = 1
265
266                 ! extrapolation pour avoir uxbar(i+1,j)
267
268                 phi_ip3 = Uxbar(i+3,j) * Hmxloc(i+3,j)
269                 phi_xgl(i+2,j) = phi_ip3 + dx * (phi_prescr - phi_ip3)/(2.5 * dx  - xpos)
270                 uxbartemp(i+2,j)   = uxbartemp(i+2,j) + phi_xgl(i+2,j) / Hmxloc(i+2,j)
271                 countx   (i+2,j)   = countx(i+2,j)+1
272                 imx_diag (i+2,j) = 1
273                 gr_line_schoof(i+2,j) = 1
274
275!afq --                 end if
276             
277              end if
278
279           end if
280
281           !*******************************************************************************
282           ! selon y (indice j)
283           if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) then 
284           !if (flot(i,j-1).and.Uybar(i,j).LT.0.) then    ! grounding line between j-1,j    CAS SUD (ecoulement vers grille SUD)
285
286
287              call  calc_xgl4schoof(-dy,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i,j-1)),real(Bsoc(i,j-1)),sealevel_2d(i,j-1),ypos,Hgl)
288              ypos_tab(i,j)=ypos
289              Hgly_tab(i,j)=Hgl
290             
291              ! afq: the back force is on the staggered grid, the GL can be either South or North to this point.
292              if (ypos .lt. -dy/2.) then
293                 bfy = back_force_y(i,j)
294              else
295                 bfy = back_force_y(i,j+1)
296              endif
297              if (gr_select.eq.1) then ! flux de Tsai
298                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr)
299              else if (gr_select.eq.2) then ! flux de Schoof
300                 ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point.
301                 if (ypos .lt. -dy/2.) then
302                    frot_coef = betamy(i,j)
303                 else
304                    frot_coef = betamy(i,j+1)
305                 endif
306                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr)
307              else
308                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
309                 stop
310              endif
311              phi_prescr=-phi_prescr ! Doit être negatif dans le cas Sud
312              phi_prescr_taby(i,j)=phi_prescr
313              if (ypos .lt. -dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag
314
315                 !  grille    ! .....x......o......x......o...G..x......O......x......o......x......o
316                 !            !            j-2           j-1            j            j+1           j+2
317                 !flux        !     in            out         G out            in
318
319!afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j-1))) then
320
321                 ! extrapolation pour avoir uybar(i,j-1) qui sera ensuite prescrit dans call rempli_L2
322
323                 phi_jm2 = Uybar(i,j-2) * Hmyloc(i,j-2)
324                 phi_ygl(i,j-1)  = phi_jm2 + dy * (phi_prescr-phi_jm2)/(2.5 * dy + ypos)
325                 uybartemp(i,j-1)   = uybartemp(i,j-1) + phi_ygl(i,j-1) / Hmyloc(i,j-1) ! attention division par 0 possible ?
326                 county   (i,j-1)   = county(i,j-1)+1
327                 imy_diag (i,j-1) = 1
328                 gr_line_schoof(i,j-1) = 1
329
330                 ! extrapolation pour avoir uybar(i,j)
331
332                 phi_jp1 = Uybar(i,j+1) * Hmyloc(i,j+1)
333                 phi_ygl(i,j) = phi_jp1 + dy * (phi_prescr - phi_jp1)/(0.5 * dy - ypos)
334                 uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j)
335                 county   (i,j)   = county(i,j)+1
336                 imy_diag (i,j) = 1
337                 gr_line_schoof(i,j) = 1
338
339
340!afq --                 end if
341             
342              else            !  GL au nord du j staggered,                          o centre, x stag
343
344                 !  grille    ! .....x......o......x......o......x...G..O......x......o......x......o
345                 !            j-2           j-1            j            j+1           j+2
346                 !flux        !     in             in           out  G         out           in
347
348!afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then
349                   
350                 ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
351
352                 phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1)
353                 phi_ygl(i,j)  = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos)
354                 uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ?
355                 county   (i,j)   = county(i,j)+1
356                 imy_diag(i,j) = 1
357                 gr_line_schoof(i,j) = 1
358
359                 ! extrapolation pour avoir uybar(i,j+1)
360
361                 phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2)
362                 phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5 * dy  - ypos)
363                 uybartemp(i,j+1)   = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1)
364                 county   (i,j+1)   = county(i,j+1)+1
365                 imy_diag (i,j+1) = 1
366                 gr_line_schoof(i,j+1) = 1
367
368!afq --                 end if
369             
370              end if
371             
372           else if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) then
373           !else if (flot(i,j+1).and.Uybar(i,j+1).GT.0.) then    ! grounding line between j,j+1    CAS NORD (ecoulement vers grille Nord)
374
375              call  calc_xgl4schoof(dy,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i,j+1)),real(Bsoc(i,j+1)),sealevel_2d(i,j+1),ypos,Hgl)
376              ypos_tab(i,j)=ypos
377              Hgly_tab(i,j)=Hgl
378             
379              ! afq: the back force is on the staggered grid, the GL can be either South or North to this point.
380              if (ypos .lt. dy/2.) then
381                 bfy = back_force_y(i,j)
382              else
383                 bfy = back_force_y(i,j+1)
384              endif
385              if (gr_select.eq.1) then ! flux de Tsai
386                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr)
387              else if (gr_select.eq.2) then ! flux de Schoof
388                 ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point.
389                 if (ypos .lt. dy/2.) then
390                    frot_coef = betamy(i,j)
391                 else
392                    frot_coef = betamy(i,j+1)
393                 endif
394                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr)
395              else
396                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
397                 stop
398              endif
399              phi_prescr_taby(i,j)=phi_prescr
400              if (ypos .lt. dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag
401
402                 !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o
403                 !            !            j-2           j-1            j            j+1           j+2
404                 !flux        !                    in           out        G  out            in
405
406!afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then
407
408                 ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
409
410                 phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1)
411                 phi_ygl(i,j)  = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos)
412                 uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ?
413                 county   (i,j)   = county(i,j)+1
414                 imy_diag (i,j) = 1
415                 gr_line_schoof(i,j) = 1
416
417                 ! extrapolation pour avoir uybar(i,j+1)
418
419                 phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2)
420                 phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5* dy - ypos)
421                 uybartemp(i,j+1)   = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1)
422                 county   (i,j+1)   = county(i,j+1)+1
423                 imy_diag (i,j+1) = 1
424                 gr_line_schoof(i,j+1) = 1
425
426!afq --                 end if
427             
428              else   !  GL au nord du j staggered,                          o centre, x stag
429
430                 !  grille    ! ......x......o......x......O......x...G..o......x......o......x.....o
431                 !            !             j-1            j            j+1           j+2          j+3
432                 !flux        !                     in           out  G        out            in
433
434!afq --                 if (.not.(ilemy(i,j+1)).and..not.(ilemy(i,j+2))) then
435
436                 ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
437
438                 phi_j = Uybar(i,j) * Hmyloc(i,j)
439                 phi_ygl(i,j+1)  = phi_j + dy * (phi_prescr-phi_j)/(0.5 * dy + ypos)
440                 uybartemp(i,j+1)   = uybartemp(i,j+1) +  phi_ygl(i,j+1) / Hmyloc(i,j+1) ! attention division par 0 possible ?
441                 county   (i,j+1)   = county(i,j+1)+1
442                 imy_diag (i,j+1) = 1
443                 gr_line_schoof(i,j+1) = 1
444
445                 ! extrapolation pour avoir uybar(i,j+1)
446
447                 phi_jp3 = Uybar(i,j+3) * Hmyloc(i,j+3)
448                 phi_ygl(i,j+2) = phi_jp3 + dy * (phi_prescr - phi_jp3)/(2.5 * dy  - ypos)
449                 uybartemp(i,j+2)   = uybartemp(i,j+2) +  phi_ygl(i,j+2) / Hmyloc(i,j+2)
450                 county   (i,j+2)   = county(i,j+2)+1
451                 imy_diag (i,j+2) = 1
452                 gr_line_schoof(i,j+2) = 1
453
454!afq --                 end if
455
456              end if
457
458           end if
459        end if
460     end do
461  end do
462
463  ! afq -- if we discard the points with multiple fluxes coming, uncomment:
464  where (countx(:,:).ge.2)
465     uxbartemp(:,:)=uxbar(:,:)*countx(:,:)
466     imx_diag(:,:) = imx_diag_in(:,:)
467  end where
468  where (county(:,:).ge.2)
469     uybartemp(:,:)=uybar(:,:)*county(:,:)
470     imy_diag(:,:) = imx_diag_in(:,:)
471  end where
472  ! afq --
473 
474  where (countx(:,:).ne.0) uxbar(:,:)= max(min(uxbartemp(:,:)/countx(:,:),V_limit),-V_limit)
475  where (county(:,:).ne.0) uybar(:,:)= max(min(uybartemp(:,:)/county(:,:),V_limit),-V_limit)
476           
477  do j= 3,ny-3
478     do i=3,nx-3
479        denom = sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2))
480        if (denom.gt.toutpetit) then
481           ! prodscal is cos theta in : U . V = u v cos theta
482           ! this number is between -1 and 1, 1 = same direction & -1 = opposite direction
483           ! we could use this as an artificial back force...
484           prodscal = (uxbartemp(i,j)*uxbar(i,j) + uybartemp(i,j)*uybar(i,j))/  &
485                (sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2)))
486        endif
487     end do
488  end do
489 
490!!$  do j=1,ny
491!!$     do i=1,nx
492!!$        write(598,*) xpos_tab(i,j)
493!!$        write(599,*) ypos_tab(i,j)
494!!$        write(600,*) Hglx_tab(i,j)
495!!$        write(601,*) Hgly_tab(i,j)
496!!$        write(602,*) Abar(i,j)
497!!$        write(603,*) phi_prescr_tabx(i,j)/Hglx_tab(i,j)
498!!$        write(604,*) phi_prescr_taby(i,j)/Hgly_tab(i,j)
499!!$     enddo
500!!$  enddo
501
502  debug_3D(:,:,66) = phi_prescr_tabx(:,:)
503  debug_3D(:,:,67) = phi_prescr_taby(:,:)
504 
505end subroutine interpol_glflux
506 
507!-------------------------------------------------------------------------
508! calcule la position sous maille de la ligne dechouage
509 subroutine calc_xgl4schoof(dy,alpha,H_0,B_0,SL_0,H_1,B_1,SL_1,xpos,Hgl)
510               
511
512  implicit none
513
514  real,intent(in)        ::    dy        !< longueur orientee de la maille
515  real,intent(in)        ::    alpha     !< coefficient de flottaison = rho/rhow
516  real,intent(in)        ::    H_0       !< epaisseur au point pose
517  real,intent(in)        ::    B_0       !< altitude socle au point pose
518  real,intent(in)        ::    SL_0      !< sea level au point pose
519  real,intent(in)        ::    H_1       !< epaisseur au point flottant
520  real,intent(in)        ::    B_1       !< altitude socle au point flottant
521  real,intent(in)        ::    SL_1      !< sea level au point flottant
522  real,intent(out)       ::    xpos      !< position de la ligne (en distance depuis le point pose)
523  real,intent(out)       ::    Hgl       !< epaisseur de glace a la grounding line
524   
525 
526
527  real                   ::    Cgl       ! variable de travail 
528
529!~   Cgl   = (-alpha * (H_1-H_0) - (B_1-B_0)) ! -(rho/rhow (H1-H0) + (B1-B0))
530
531!~   if (abs(Cgl).gt.1.e-5) then   
532!~      xpos    = (B_0 + alpha * H_0 - SL) / Cgl
533!~   else
534!~      xpos    = 1. - 1./abs(dy)   ! verifier
535!~   end if
536 
537!~   if (xpos.LT.0..OR.xpos.GT.1.) then
538!~              print*,'calc_xgl4schoof : xpos value error i,j=',i,j
539!~              print*, 'xpos,Cgl', xpos,Cgl
540!~              print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL
541!~              stop
542!~      endif
543               
544  Cgl   = (B_1-B_0) + alpha * (H_1-H_0) - (SL_1-SL_0)
545
546  if (abs(Cgl).gt.1.e-5) then   
547     xpos    = (SL_0 - (B_0 + alpha * H_0)) / Cgl 
548  else
549     xpos    = 0.5  ! verifier
550  end if
551 
552!  print*
553!  print*, 'xpos', xpos
554!       print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL
555!       print*,'B_1, H_1', B_1, H_1
556
557       
558  if (xpos.LT.0..OR.xpos.GT.1.) then
559     print*,'calc_xgl4schoof : xpos value error i,j=',i,j
560     print*, 'xpos,Cgl', xpos,Cgl
561     print*,'B_0, alpha, H_0, SL_0', B_0, alpha, H_0, SL_0
562     print*,'archim=',B_1+H_1*alpha - SL_1
563     !stop
564     xpos = min(max(0.,xpos),1.)
565  endif
566       
567  xpos = xpos * dy
568 
569  Hgl= H_0 + (H_1 - H_0) * xpos/dy
570
571end subroutine calc_xgl4schoof
572
573
574!--------------------------------------------------------------------
575! Schoof_flux : calcule le coefficient et le flux de Schoof
576!--------------------------------------------------------------------
577
578! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement)
579
580subroutine flux_Schoof4Schoof (Hgl,A_mean,C_frot,alpha,n_glen,m_weert,back_force_coef,phi_schoof)
581
582 
583implicit none
584
585real,intent(in) :: Hgl
586real,intent(in) :: A_mean
587real,intent(in) :: C_frot
588real,intent(in) :: alpha
589real,intent(in) :: n_glen
590real,intent(in) :: m_weert
591real,intent(in) :: back_force_coef
592real,intent(out) :: phi_schoof
593
594! afq: in standard GRISLI (19/05/17), tob= -beta ub => m_weert = 1 and C_frot = beta
595
596!phi_schoof = (((A_mean * rog**(n_glen+1) * (1 - alpha)**n_glen) / (4**n_glen *C_frot)) **(1./(1.+1./m_weert)))*back_force_coef**(n_glen/(1+1./m_weert))
597!phi_schoof = phi_schoof * Hgl**((n_glen+3.+1/m_weert)/(1.+1./m_weert))
598phi_schoof = (((A_mean * rog**(n_glen+1.) * (1. - alpha)**n_glen) / (4.**n_glen*C_frot)) **(1./(1.+inv_mweert)))*back_force_coef**(n_glen/(1.+inv_mweert))
599phi_schoof = phi_schoof * Hgl**((n_glen+3.+inv_mweert)/(1.+inv_mweert))
600
601
602
603end subroutine flux_Schoof4Schoof
604
605!--------------------------------------------------------------------
606! Tsai_flux : calcule le coefficient et le flux de Schoof
607!--------------------------------------------------------------------
608
609! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement)
610
611subroutine flux_Tsai4Schoof (Hgl,A_mean,f_frot,alpha,n_glen,back_force_coef,phi_Tsai)
612
613implicit none
614
615real,intent(in) :: Hgl
616real,intent(in) :: A_mean
617real,intent(in) :: f_frot
618real,intent(in) :: alpha
619real,intent(in) :: n_glen
620real,intent(in) :: back_force_coef
621real,intent(out) :: phi_Tsai
622
623real, parameter :: Q0 = 0.61
624
625
626phi_Tsai = Q0 * ((8. * A_mean * rog**(n_glen) * (1. - alpha)**(n_glen-1.)) / (4.**n_glen *f_frot))*back_force_coef**(n_glen-1.)
627phi_Tsai = phi_Tsai * Hgl**(n_glen+2.)
628
629end subroutine flux_Tsai4Schoof
630
631end module furst_schoof_mod
Note: See TracBrowser for help on using the repository browser.