source: trunk/SOURCES/furst_schoof_mod.f90 @ 175

Last change on this file since 175 was 175, checked in by dumas, 6 years ago

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