source: trunk/SOURCES/furst_schoof_mod.f90 @ 176

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

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