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