- Timestamp:
- 02/09/24 17:08:19 (4 months ago)
- Location:
- branches/GRISLIv3/SOURCES/New-remplimat
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/GRISLIv3/SOURCES/New-remplimat/diagno-L2_mod.f90
r481 r482 16 16 real, dimension(nx,ny) :: uyb1ramollo 17 17 real, dimension(nx,ny) :: pvi_keep 18 19 !cdc transfere dans module3d pour compatibilite avec furst_schoof_mod20 !cdc integer, dimension(nx,ny) :: imx_diag21 !cdc integer, dimension(nx,ny) :: imy_diag22 18 23 19 integer :: nneigh=max(int(400000./dx),1) !nb. points around the grline kept for back force comp. in case we reduce the extent … … 151 147 !$OMP END WORKSHARE 152 148 !$OMP END PARALLEL 153 ! avec couplage thermomecanique154 ! write(166,*) ' apres call calc_pvi',itour_pvi155 149 156 150 else … … 158 152 end if 159 153 160 call imx_imy_nx_ny ! pour rempli_L2 : calcule les masques imx et imy qui154 call imx_imy_nx_ny(imx_diag,imy_diag,flgzmx,flgzmy) ! pour rempli_L2 : calcule les masques imx et imy qui 161 155 162 156 diagno_grline=sum(gr_line(2:(nx-1),2:(ny-1))) 163 157 if ((Schoof.ge.1).and.(diagno_grline.gt.0)) then 164 call imx_imy_nx_ny_reduce(1 ) !afq: provides imx_reduce and imy_reduce158 call imx_imy_nx_ny_reduce(1,flot,imx_diag,imy_diag,gr_line) !afq: provides imx_reduce and imy_reduce 165 159 endif 166 167 !cdc debug Schoof !!!!!!!!!!!!168 !~ do j=1,ny169 !~ do i=1,nx170 !~ write(578,*) uxbar(i,j)171 !~ write(579,*) uybar(i,j)172 !~ enddo173 !~ enddo174 175 !if (Schoof.eq.1.and.nt.GT.15000) then ! flux grounding line Schoof176 ! afq -- below: if (Schoof.eq.1) then ! flux grounding line Schoof177 ! afq -- below: call interpol_glflux ! calcul flux GL + interpolation sur voisins178 ! afq -- below: endif179 180 !~ do j=1,ny181 !~ do i=1,nx182 !~ write(588,*) uxbar(i,j)183 !~ write(589,*) uybar(i,j)184 !~ enddo185 !~ enddo186 !~ print*,'ecriteure termineee !!!!!!'187 !~ read(*,*)188 160 189 161 ! donnent les cas de conditions aux limites … … 202 174 nyd1=1 203 175 nyd2=ny 204 205 !call rempli_L2(1,nx,1,ny,uxbar,uybar,uxb1,uyb1,imx_diag,imy_diag,ifail_diagno)206 !nxd1=15207 !nxd2=19208 !nyd1=30209 !nyd2=34210 211 !nxd1=35212 !nxd2=60213 !nyd1=35214 !nyd2=60215 216 176 217 !if (Schoof.eq.1.and.nt.GT.15000) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo218 177 if ((Schoof.ge.1).and.(diagno_grline.gt.0)) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo 219 178 … … 224 183 pvi_keep(:,:)=pvi(:,:) 225 184 where (flot(:,:).and.H(:,:).GT.2.) 226 ! pvi(:,:)=1.e5227 185 pvi(:,:)=pvimin 228 186 endwhere … … 250 208 251 209 if (ifail_diagno_ramollo.gt.0) then 252 ! write(6,*) ' Probleme resolution systeme L2. ramollo ifail=',ifail_diagno_ramollo253 ! STOP254 210 write(*,*) ' Probleme resolution systeme L2. ramollo ifail=',ifail_diagno_ramollo 255 211 write(*,*) ' ... we go on anyway!' 256 212 endif 257 !~ do j=1,ny258 !~ do i=1,nx259 !~ if (sqrt(uxb1(i,j)**2+ uyb1(i,j)**2).gt.0..and..not.flot(i,j)) then260 !~ write(1034,*) sqrt(uxb1(i,j)**2+ uyb1(i,j)**2) / sqrt(uxb1ramollo(i,j)**2 + uyb1ramollo(i,j)**2)261 !~ else262 !~ write(1034,*) 1.263 !~ endif264 !~ enddo265 !~ enddo266 267 !~ print*,'apres calcul rempli_L2'268 !~ read(*,*)269 213 270 214 call interpol_glflux ! calcul flux GL + interpolation sur voisins … … 503 447 debug_3D(:,:,27)=pvi(:,:) 504 448 !$OMP END WORKSHARE 505 ! attention run 35506 !--------------------507 !!$if (time.gt.10.) then508 !!$ where (flot(:,:))509 !!$ pvi(:,:)=pvimin510 !!$ end where511 !!$end if512 449 513 450 ! calcul de la viscosite integree au milieu des mailles (pvm) … … 528 465 !------------------------------------------------------------------ 529 466 530 subroutine imx_imy_nx_ny 531 532 use module3d_phy, only: imx_diag,imy_diag,flgzmx,flgzmy 533 467 subroutine imx_imy_nx_ny(imx_diagl,imy_diagl,flgzmxl,flgzmyl) 468 534 469 implicit none 470 471 logical,dimension(:,:), intent(in) :: flgzmxl,flgzmyl 472 integer,dimension(:,:), intent(inout) :: imx_diagl,imy_diagl 535 473 536 474 ! definition des masques … … 555 493 !$OMP PARALLEL 556 494 !$OMP WORKSHARE 557 imx_diag (:,:)=0558 imy_diag (:,:)=0495 imx_diagl(:,:)=0 496 imy_diagl(:,:)=0 559 497 560 498 ! a l'interieur du domaine 561 499 !------------------------- 562 500 563 where (flgzmx (:,:))564 imx_diag (:,:)=2 ! shelf ou stream501 where (flgzmxl(:,:)) 502 imx_diagl(:,:)=2 ! shelf ou stream 565 503 elsewhere 566 imx_diag (:,:)=1 ! vitesse imposee504 imx_diagl(:,:)=1 ! vitesse imposee 567 505 end where 568 506 569 where (flgzmy (:,:))570 imy_diag (:,:)=2 ! shelf ou stream507 where (flgzmyl(:,:)) 508 imy_diagl(:,:)=2 ! shelf ou stream 571 509 elsewhere 572 imy_diag (:,:)=1 ! vitesse imposee510 imy_diagl(:,:)=1 ! vitesse imposee 573 511 end where 574 512 575 513 ! bord sud 576 imx_diag (:,1)=-1577 imy_diag (:,2)=-1514 imx_diagl(:,1)=-1 515 imy_diagl(:,2)=-1 578 516 579 517 ! bord nord 580 imx_diag (:,ny)=-3581 imy_diag (:,ny)=-3518 imx_diagl(:,ny)=-3 519 imy_diagl(:,ny)=-3 582 520 583 521 ! bord Est 584 imx_diag (1,:)=0 ! hors domaine a cause des mailles alternees585 imx_diag (2,:)=-4586 imy_diag (1,:)=-4522 imx_diagl(1,:)=0 ! hors domaine a cause des mailles alternees 523 imx_diagl(2,:)=-4 524 imy_diagl(1,:)=-4 587 525 588 526 ! bord West 589 imx_diag (nx,:)=-2590 imy_diag (nx,:)=-2527 imx_diagl(nx,:)=-2 528 imy_diagl(nx,:)=-2 591 529 592 530 ! Coins 593 imx_diag (2,1)=-41 ! SW594 imy_diag (1,2)=-41595 596 imx_diag (nx,1)=-12 ! SE597 imy_diag (nx,2)=-12598 599 imx_diag (nx,ny)=-23 ! NE600 imy_diag (nx,ny)=-23601 602 imx_diag (2,ny)=-34 ! NW603 imy_diag (1,ny)=-34531 imx_diagl(2,1)=-41 ! SW 532 imy_diagl(1,2)=-41 533 534 imx_diagl(nx,1)=-12 ! SE 535 imy_diagl(nx,2)=-12 536 537 imx_diagl(nx,ny)=-23 ! NE 538 imy_diagl(nx,ny)=-23 539 540 imx_diagl(2,ny)=-34 ! NW 541 imy_diagl(1,ny)=-34 604 542 605 543 ! hors domaine 606 imx_diag (1,:)=0 ! hors domaine a cause des mailles alternees607 imy_diag (:,1)=0 ! hors domaine a cause des mailles alternees544 imx_diagl(1,:)=0 ! hors domaine a cause des mailles alternees 545 imy_diagl(:,1)=0 ! hors domaine a cause des mailles alternees 608 546 !$OMP END WORKSHARE 609 547 !$OMP END PARALLEL … … 612 550 613 551 614 subroutine imx_imy_nx_ny_reduce(choix) 615 616 use module3d_phy, only: flot,imx_diag,imy_diag,gr_line 617 552 subroutine imx_imy_nx_ny_reduce(choix,flot,imx_diag,imy_diag,gr_line) 553 618 554 implicit none 619 555 … … 621 557 ! We simply compute the velocities around (nneigh) the grounding line 622 558 623 integer, intent(in) :: choix 559 integer, intent(in) :: choix 560 logical, dimension(:,:), intent(in) :: flot 561 integer, dimension(:,:), intent(in) :: imx_diag,imy_diag,gr_line 562 624 563 625 564 integer i,j,nvx,nvy … … 717 656 end subroutine mismip_boundary_cond 718 657 719 !end module diagno_L2_mod720 658 end module diagno_mod -
branches/GRISLIv3/SOURCES/New-remplimat/eq_ellipt_sgbsv_mod-0.2.f90
r478 r482 57 57 58 58 !---------------------------- 59 subroutine redim_matrice 59 subroutine redim_matrice(mmatl) 60 60 61 61 ! allocation dynamique de mmat … … 63 63 implicit none 64 64 65 real,dimension(:,:), allocatable, intent(inout) :: mmatl 66 65 67 integer :: err 66 68 67 if (allocated(mmat )) then68 deallocate(mmat ,stat=err)69 if (allocated(mmatl)) then 70 deallocate(mmatl,stat=err) 69 71 if (err/=0) then 70 72 print *,"Erreur à la desallocation de mmat",err … … 74 76 75 77 76 allocate(mmat (ldb,nblig),stat=err)78 allocate(mmatl(ldb,nblig),stat=err) 77 79 if (err/=0) then 78 80 print *,"erreur a l'allocation du tableau mmat",err … … 80 82 stop 81 83 end if 82 83 84 !write(6,*) "allocation mmat ok, ldb,nblig",ldb,nblig85 84 86 85 end subroutine redim_matrice … … 364 363 ! si la matrice mmat change de taille on la re allocate. 365 364 366 if ((ldb.ne.ldb_old).and.(nblig.ne.nblig_old)) call redim_matrice 365 if ((ldb.ne.ldb_old).and.(nblig.ne.nblig_old)) call redim_matrice(mmat) 367 366 368 367 -
branches/GRISLIv3/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90
r481 r482 131 131 !$OMP END PARALLEL 132 132 133 ! calcul de Hmx_oppos et Hmy_oppos dans le cas de fronts134 ! ce calcul pourrait être dans diagno et faire passer hmx_oppos par module135 136 ! 4 avril 2012 : Finalement il me semble que c'est faux, il suffit de garder le Hmx habituel.137 ! cat138 139 !!$do j=ny1,ny2140 !!$ do i=nx1,nx2141 !!$142 !!$ il=max(1,i-1)143 !!$ jl=max(1,j-1)144 !!$145 !!$146 !!$ if ((H(i,j).gt.1.).and.(H(il,j).le.1.)) then ! Bord West : on garde H(i,j)147 !!$ Hmx_oppos(i,j)=H(i,j)148 !!$149 !!$ else if ((H(i,j).le.1.).and.(H(il,j).gt.1.)) then ! Bord Est : on garde H(il,j)150 !!$ Hmx_oppos(i,j)=H(il,j)151 !!$152 !!$ end if153 !!$154 !!$ if ((H(i,j).gt.1.).and.(H(i,jl).le.1.)) then ! Bord Sud : on garde H(i,j)155 !!$ Hmy_oppos(i,j)=H(i,j)156 !!$157 !!$ else if ((H(i,j).le.1.).and.(H(i,jl).gt.1.)) then ! Bord Nord on garde H(i,jl)158 !!$ Hmy_oppos(i,j)=H(i,jl)159 !!$160 !!$ end if161 !!$162 !!$163 !!$ end do164 !!$ end do165 166 167 ! mailles mineures en dehors du domaine168 !imy(:,ny1) = 0169 !imx(nx1,:) = 0170 171 133 172 134 ! limite la valeur de frotmx a betamax … … 175 137 call limit_frotm 176 138 177 178 139 ! remplissage des Tuij, .... 179 140 !---------------------------- 180 141 call rempli_Tuij 181 182 ! debug : impression des Tuij en un point183 !write(166,*)184 !write(166,*) 'impression des Tuij'185 186 !j = 12187 188 !do i = 244,245189 190 !!$do i = 244,246191 !!$ write(166,801) i,imx(i,j),frotmx(i,j),hmx_oppos(i,j),sdx(i,j)192 !!$ do k= 1,-1,-1193 !!$ write(166,800) k,Tu(i,j,-1:1,k),opposx(i,j)194 !!$ end do195 !!$ write(166,*)196 !!$ do k= 1,-1,-1197 !!$ write(166,800) k,Tv(i,j,-1:1,k)198 !!$ end do199 !!$ write(166,*)200 !!$ write(166,*) 'impression des Svij'201 !!$ write(166,801) i,imy(i,j),frotmy(i,j),hmy_oppos(i,j),sdy(i,j)202 !!$ do k= 1,-1,-1203 !!$ write(166,800) k,Su(i,j,-1:1,k),opposy(i,j)204 !!$ end do205 !!$ write(166,*)206 !!$ do k= 1,-1,-1207 !!$ write(166,800) k,Sv(i,j,-1:1,k)208 !!$ end do209 !!$210 !!$211 !!$212 !!$800 format(' k=',(i0,1x),': ',4(es12.4,1x) )213 !!$801 format(' i,im,frotm,hm,sd ',2(i0.1,x),3(es12.4,1x))214 !!$ write(166,*)215 !!$end do216 !!$write(166,*)217 218 219 ! debug : differents termes de l equation MacAyeal220 !---------------------------------------------------221 222 223 !!$do j=ny1+1,ny2-1224 !!$ do i=nx1+1,nx2-1225 !!$226 !!$ if (gzmx(i,j)) then227 !!$ do jl=-1,1228 !!$ do il=-1,1229 !!$ debug_3D(i,j,82) =debug_3D(i,j,82)+(Tu(i,j,il,jl)*Vcol_x(i+il,j+jl) &230 !!$ +Tv(i,j,il,jl)*Vcol_y(i+il,j+jl))231 !!$232 !!$ end do233 !!$ end do234 !!$! debug_3D(i,j,83)=frotmx(i,j)*dx2*Vcol_x(i,j)235 !!$ end if236 !!$ end do237 !!$end do238 239 240 !!$do j=ny1,ny2241 !!$ do i=nx1,nx2242 !!$ debug_3D(i,j,84) = imx(i,j)243 !!$ end do244 !!$end do245 246 ! call rempli_Tuij ! pour reutiliser le nouveau beta247 248 142 249 143 ! Pour determiner les colonnes Mu, Mv, Nu, Nv … … 268 162 if(abs(scal).lt.1.e-30) then 269 163 write(num_tracebug,*)' i,j,imx,pvi',i,j,imx(i,j),pvi(i,j) 270 ! write(num_tracebug,*) Tu(i,j,:,:)271 164 end if 272 165 … … 275 168 Tv(i,j,:,:)=Tv(i,j,:,:)/scal ! Tu(:,:,0,0) est la diagonale 276 169 opposx(i,j)=opposx(i,j)/scal 277 278 ! if ((i.eq.200).and.(j.eq.251)) then279 ! write(num_tracebug,*)' i,j,imx',i,j,imx(i,j)280 ! write(num_tracebug,*)' pvi,scal',pvi(i,j),scal281 ! write(num_tracebug,*) Tu(i,j,:,-1)282 ! write(num_tracebug,*) Tu(i,j,:,0)283 ! write(num_tracebug,*) Tu(i,j,:,1)284 ! write(num_tracebug,*)' opposx',opposx(i,j)285 ! end if286 287 288 170 289 171 end if … … 293 175 if(abs(scal).lt.1.e-30) then 294 176 write(num_tracebug,*)' i,j,Svij,imy,pvi',i,j,imy(i,j),pvi(i,j) 295 ! write(num_tracebug,*) Sv(i,j,:,:)296 177 end if 297 178 … … 339 220 !$OMP END PARALLEL 340 221 341 ! pour appeler la routine de resolution du systeme lineaire342 343 !!$i = 244344 !!$do j=11,11,-1345 !!$346 !!$write(166,*) 'j=',j347 !!$ write(166,803) i,uxprec(i,j),i+1,uxprec(i+1,j),i+2,uxprec(i+2,j)348 !!$ write(166,802) i,uxnew(i,j),i+1,uxnew(i+1,j),i+2,uxnew(i+2,j)349 !!$802 format(3('uxnew ',' i=',(i0,1x,es12.4,1x) ))350 !!$803 format(3('uxprec',' i=',(i0,1x,es12.4,1x) ))351 !!$write(166,*)352 353 !!$end do354 355 ! limitation des nouvelles vitesses356 357 !!$uxnew(:,:)=max(-3900.,uxnew(:,:))358 !!$uxnew(:,:)=min( 3900.,uxnew(:,:))359 !!$uynew(:,:)=max(-3900.,uynew(:,:))360 !!$uynew(:,:)=min( 3900.,uynew(:,:))361 !!$362 363 364 365 366 222 888 return 367 368 223 369 224 … … 695 550 real, intent(out) :: opposx 696 551 697 698 ! Tu(i,j,0,0) = 4.*pvi(i-1,j)699 ! Tu(i,j,1,0) = -4.*pvi(i-1,j)700 701 552 Tu(0,0) = 4.*pvi ! quand le bord est epais, il faut prendre 702 553 Tu(1,0) = -4.*pvi ! le noeud du shelf pas celui en H=1 703 ! s'il est fin, pas de difference 704 705 !!$ Tv(i-1,j,-1,0) = 2.*pvi(i-1,j) ! ???? 706 !!$ Tv(i-1,j,-1,1) = -2.*pvi(i-1,j) 554 ! s'il est fin, pas de difference 707 555 708 556 Tv(-1,0) = 2.*pvi_imoinsun … … 1227 1075 !$OMP END WORKSHARE 1228 1076 !$OMP END PARALLEL 1229 1230 ! sortie Netcdf pour verifier ok_umat1231 !~ where (ok_umat(:,:))1232 !~ debug_3D(:,:,39)=11233 !~ elsewhere1234 !~ debug_3D(:,:,39)=01235 !~ end where1236 !~1237 !~ where (ok_vmat(:,:))1238 !~ debug_3D(:,:,40)=11239 !~ elsewhere1240 !~ debug_3D(:,:,40)=01241 !~ end where1242 !~1243 !~ ! sortie Netcdf pour verifier ghost1244 !~ where (ghost_x(:,:))1245 !~ debug_3D(:,:,41)=11246 !~ elsewhere1247 !~ debug_3D(:,:,41)=01248 !~1249 !~ end where1250 !~1251 !~ where (ghost_y(:,:))1252 !~ debug_3D(:,:,42)=11253 !~ elsewhere1254 !~ debug_3D(:,:,42)=01255 !~ end where1256 1077 1257 1078 if (itracebug.eq.1) call tracebug(' Sortie Okmat') … … 1374 1195 do il=-1,1 1375 1196 betamx(i,j) = betamx(i,j)+(Tu(i,j,il,jl)*uxgiven(i+il,j+jl) & 1376 +Tv(i,j,il,jl)*uygiven(i+il,j+jl)) 1377 1378 ! debug_3D(i,j,82)=debug_3D(i,j,82) + (Tu(i,j,il,jl)*uxgiven(i+il,j+jl) & 1379 ! +Tv(i,j,il,jl)*uygiven(i+il,j+jl)) 1380 1197 +Tv(i,j,il,jl)*uygiven(i+il,j+jl)) 1381 1198 end do 1382 1199 end do 1383 1384 ! debug_3D(i,j,82)=debug_3D(i,j,82)/dx/dx ! longitudinal stress computed from velocities1385 1386 1200 1387 1201 ! fonctionne même quand frotmx n'est pas nul … … 1408 1222 1409 1223 betamy(i,j)=-opposy(i,j) 1410 ! debug_3D(i,j,36) = opposy(i,j)/dx/dx ! driving stress ou pression hydrostatique (Pa)1411 1224 do jl=-1,1 1412 1225 do il=-1,1 … … 1415 1228 +Sv(i,j,il,jl)*uygiven(i+il,j+jl)) 1416 1229 1417 ! debug_3D(i,j,83)=debug_3D(i,j,83) + (Su(i,j,il,jl)*uxgiven(i+il,j+jl) &1418 ! +Sv(i,j,il,jl)*uygiven(i+il,j+jl))1419 1420 1421 1230 end do 1422 1231 end do 1423 1424 ! debug_3D(i,j,83)=debug_3D(i,j,83) /dx/dx ! longitudinal stress computed from velocities1425 1232 1426 1233 ! fonctionne même quand frotmx n'est pas nul … … 1479 1286 !$OMP END PARALLEL 1480 1287 1481 1482 ! les noeuds negatifs veulent dire qu'il faudrait ajouter un moteur pour aller aussi vite1483 ! que les vitesses de bilan.1484 ! il faudrait repartir plutot que mettre a 01485 1486 ! a mettre ailleurs1487 ! betamx(:,:)=max(betamx(:,:),0.)1488 ! betamy(:,:)=max(betamy(:,:),0.)1489 ! beta_centre(:,:)=max(beta_centre(:,:),0.)1490 1491 ! iter_beta=21492 1493 1494 1288 if (itracebug.eq.1) call tracebug(' Fin calc_beta') 1495 1289 end subroutine calc_beta
Note: See TracChangeset
for help on using the changeset viewer.