source: trunk/SOURCES/furst_schoof_mod.f90 @ 170

Last change on this file since 170 was 110, checked in by aquiquet, 7 years ago

Implementation of Schoof4Schoof parametrisation instead of Tsai4Schoof. Detection of floating points stuck inside grounded points (marais). Islands and marais are discarded for Schoof computation. flottab reverted to determin_tache and front (icebergs have to be put in calving for mass conservation).

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