Changeset 192 for trunk/SOURCES/New-remplimat/diagno-L2_mod.f90

Timestamp:

05/29/18 17:49:59 (6 years ago)

Author:

aquiquet

Message:

Model performance improved with a reduced domain extension for elliptic equation resolution for buttressing estimation

File:

• trunk/SOURCES/New-remplimat/diagno-L2_mod.f90 (modified) (7 diffs)

trunk/SOURCES/New-remplimat/diagno-L2_mod.f90

@@ -20 +20 @@
 !cdc  integer, dimension(nx,ny) :: imx_diag
 !cdc  integer, dimension(nx,ny) :: imy_diag
+
+integer :: nneigh=max(int(400000./dx),1) !nb. points around the grline kept for back force comp. in case we reduce the extent
+integer, dimension(nx,ny) :: imx_reduce !afq -- to limit the number of points where vel is computed for backforce
+integer, dimension(nx,ny) :: imy_reduce !afq -- to limit the number of points where vel is computed for backforce
 
 integer :: nxd1,nxd2     ! domaine selon x Dans l'appel rempli_L2
@@ -92 +96 @@
 
 
+  integer :: diagno_grline
+
+  
   if (itracebug.eq.1)  call tracebug(' Entree dans diagnoshelf')
 
@@ -130 +137 @@
 
   call imx_imy_nx_ny         ! pour rempli_L2 : calcule les masques imx et imy qui 
+
+  diagno_grline=sum(gr_line(2:(nx-1),2:(ny-1)))
+  if ((Schoof.ge.1).and.(diagno_grline.gt.0)) then
+       call imx_imy_nx_ny_reduce(1) !afq: provides imx_reduce and imy_reduce
+  endif
 
 !cdc debug Schoof !!!!!!!!!!!!  
@@ -180 +192 @@
   !nyd2=60
 
-  call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), & 
-       uxb1(nxd1:nxd2,nyd1:nyd2),uyb1(nxd1:nxd2,nyd1:nyd2),  &
-       imx_diag(nxd1:nxd2,nyd1:nyd2),imy_diag(nxd1:nxd2,nyd1:nyd2),ifail_diagno)
        
   !if (Schoof.eq.1.and.nt.GT.15000) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo
-  if (Schoof.eq.1) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo
+  if ((Schoof.ge.1).and.(diagno_grline.gt.0)) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo
+
+    call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), & 
+          uxb1(nxd1:nxd2,nyd1:nyd2),uyb1(nxd1:nxd2,nyd1:nyd2),  &
+          imx_reduce(nxd1:nxd2,nyd1:nyd2),imy_reduce(nxd1:nxd2,nyd1:nyd2),ifail_diagno)
+
     pvi_keep(:,:)=pvi(:,:)       
     where (flot(:,:).and.H(:,:).GT.2.)
@@ -194 +208 @@
     call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), & 
        uxb1ramollo(nxd1:nxd2,nyd1:nyd2),uyb1ramollo(nxd1:nxd2,nyd1:nyd2),  &
-       imx_diag(nxd1:nxd2,nyd1:nyd2),imy_diag(nxd1:nxd2,nyd1:nyd2),ifail_diagno_ramollo)
+       imx_reduce(nxd1:nxd2,nyd1:nyd2),imy_reduce(nxd1:nxd2,nyd1:nyd2),ifail_diagno_ramollo)
        
     pvi(:,:)=pvi_keep(:,:) 
@@ -234 +248 @@
     call interpol_glflux ! calcul flux GL + interpolation sur voisins
     
-    call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), & 
+  endif
+  
+  call rempli_L2(nxd1,nxd2,nyd1,nyd2,uxbar(nxd1:nxd2,nyd1:nyd2),uybar(nxd1:nxd2,nyd1:nyd2), & 
        uxb1(nxd1:nxd2,nyd1:nyd2),uyb1(nxd1:nxd2,nyd1:nyd2),  &
        imx_diag(nxd1:nxd2,nyd1:nyd2),imy_diag(nxd1:nxd2,nyd1:nyd2),ifail_diagno)
-  
-  endif
   
 
@@ -559 +573 @@
  
 end subroutine imx_imy_nx_ny
+
+
+subroutine imx_imy_nx_ny_reduce(choix)
+
+!afq -- For the backforce computation we do not need to compute the velocities everywhere
+!       We simply compute the velocities around (nneigh) the grounding line
+
+  integer, intent(in) :: choix
+
+  integer i,j,nvx,nvy
+  integer bordouest,bordest,bordsud,bordnord
+
+  if (choix.eq.1) then ! we reduce the extent of velocity computation
+
+  imx_reduce(:,:) = 1
+  imy_reduce(:,:) = 1
+
+  where (flot(:,:)) 
+     imx_reduce(:,:) = imx_diag(:,:)
+     imy_reduce(:,:) = imy_diag(:,:)
+  endwhere
+
+  do j=1,ny
+     do i=1,nx
+        if (gr_line(i,j).eq.1) then
+           bordouest=max(i-nneigh,1)
+           bordest=min(i+nneigh,nx)
+           bordsud=max(j-nneigh,1)
+           bordnord=min(j+nneigh,ny)
+           do nvx=bordouest,bordest
+              do nvy=bordsud,bordnord
+                 imx_reduce(nvx,nvy)=imx_diag(nvx,nvy)
+                 imy_reduce(nvx,nvy)=imy_diag(nvx,nvy)
+              enddo
+           enddo
+        endif
+     enddo
+  enddo
+
+! -- we need to keep imx_diag for the domain edges for land terminating geometries
+  imx_reduce(:,1) = imx_diag(:,1)
+  imx_reduce(:,2) = imx_diag(:,2)
+  imx_reduce(:,ny) = imx_diag(:,ny)
+  imx_reduce(1,:) = imx_diag(1,:)
+  imx_reduce(2,:) = imx_diag(2,:)
+  imx_reduce(nx,:) = imx_diag(nx,:)
+
+  imy_reduce(:,1) = imy_diag(:,1)
+  imy_reduce(:,2) = imy_diag(:,2)
+  imy_reduce(:,ny) = imy_diag(:,ny)
+  imy_reduce(1,:) = imy_diag(1,:)
+  imy_reduce(2,:) = imy_diag(2,:)
+  imy_reduce(nx,:) = imy_diag(nx,:)
+
+
+  else ! we do not reduce the extent
+     imx_reduce(:,:)=imx_diag(:,:)
+     imy_reduce(:,:)=imy_diag(:,:)
+  endif
+
+end subroutine imx_imy_nx_ny_reduce
+
+
 !___________________________________________________________________________
 ! pour imposer les conditions de mismip sur les bords du fleuve