From 1abaff07fabc582e2ddaaf4977e337c85e5e9dfa Mon Sep 17 00:00:00 2001
From: Ted Mansell <ted.mansell@noaa.gov>
Date: Sun, 29 Oct 2023 15:42:36 -0500
Subject: [PATCH]  module_mp_nssl_2mom.F90: fix bug when nz > 128 where
 sedimentation did not work for k > 128

---
 physics/module_mp_nssl_2mom.F90 | 152 ++++++++++++++++++++++++--------
 1 file changed, 115 insertions(+), 37 deletions(-)

diff --git a/physics/module_mp_nssl_2mom.F90 b/physics/module_mp_nssl_2mom.F90
index a88ffe053..ad90ec81f 100644
--- a/physics/module_mp_nssl_2mom.F90
+++ b/physics/module_mp_nssl_2mom.F90
@@ -4274,7 +4274,7 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, &
 !      real gz(-nor+ng1:nz+nor),z1d(-nor+ng1:nz+nor,4)
       real dtp
       real xfall(nx,ny,na)  ! array for stuff landing on the ground
-      real xfall0(nx,ny)    ! dummy array
+!      real xfall0(nx,ny)    ! dummy array
       integer infdo
       integer jslab ! which line of xfall to use
             
@@ -4282,47 +4282,81 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, &
       real tmp, vtmax, dtptmp, dtfrac
       real, parameter :: dz = 200.
 
-      real :: xvt(nz+1,nx,3,lc:lhab) ! (nx,nz,2,lc:lhab) ! 1=mass-weighted, 2=number-weighted
-      real :: tmpn(-nor+ng1:nx+nor,-nor+ng1:ny+nor,-norz+ng1:nz+norz)
-      real :: tmpn2(-nor+ng1:nx+nor,-nor+ng1:ny+nor,-norz+ng1:nz+norz)
-      real :: z(-nor+ng1:nx+nor,-norz+ng1:nz+norz,lr:lhab)
-      real :: db1(nx,nz+1),dtz1(nz+1,nx,0:1),dz2dinv(nz+1,nx),db1inv(nx,nz+1)
+!      real :: xvt(nz+1,nx,3,lc:lhab) ! (nx,nz,2,lc:lhab) ! 1=mass-weighted, 2=number-weighted
+!      real :: tmpn(-nor+ng1:nx+nor,-nor+ng1:ny+nor,-norz+ng1:nz+norz)
+!      real :: tmpn2(-nor+ng1:nx+nor,-nor+ng1:ny+nor,-norz+ng1:nz+norz)
+!      real :: z(-nor+ng1:nx+nor,-norz+ng1:nz+norz,lr:lhab)
+!      real :: db1(nx,nz+1),dtz1(nz+1,nx,0:1),dz2dinv(nz+1,nx),db1inv(nx,nz+1)
       
-      real :: rhovtzx(nz,nx)
+!      real :: rhovtzx(nz,nx)
+
+      real, allocatable :: db1(:,:), dtz1(:,:,:),dz2dinv(:,:),db1inv(:,:) ! db1(nx,nz+1),dtz1(nz+1,nx,0:1),dz2dinv(nz+1,nx),db1inv(nx,nz+1)
+      real, allocatable :: rhovtzx(:,:)
+      real, allocatable :: xfall0(:,:), xvt(:,:,:,:),tmpn(:,:,:),tmpn2(:,:,:),z(:,:,:)
       
       double precision :: timesed1,timesed2,timesed3, zmaxsed,timesetvt,dummy
       double precision :: dt1,dt2,dt3,dt4
 
-      integer,parameter :: ngs = 128 
+      integer :: ngs ! = 512
       integer :: ngscnt,mgs,ipconc0
       
-      real ::  qx(ngs,lv:lhab) 
-      real ::  qxw(ngs,ls:lhab) 
-      real ::  cx(ngs,lc:lhab) 
-      real ::  xv(ngs,lc:lhab) 
-      real ::  vtxbar(ngs,lc:lhab,3) 
-      real ::  xmas(ngs,lc:lhab) 
-      real ::  xdn(ngs,lc:lhab) 
-      real ::  xdia(ngs,lc:lhab,3) 
-      real ::  vx(ngs,li:lhab) 
-      real ::  alpha(ngs,lc:lhab) 
-      real ::  zx(ngs,lr:lhab) 
-      logical :: hasmass(nx,lc+1:lhab)
-
-      integer igs(ngs),kgs(ngs)
-      
-      real rho0(ngs),temcg(ngs)
-
-      real temg(ngs)
-      
-      real rhovt(ngs)
-      
-      real cwnc(ngs),cinc(ngs)
-      real fadvisc(ngs),cwdia(ngs),cipmas(ngs)
-      
-      real cimasn,cimasx,cnina(ngs),cimas(ngs)
-      
-      real cnostmp(ngs)
+!     real ::  qx(ngs,lv:lhab) 
+!     real ::  qxw(ngs,ls:lhab) 
+!     real ::  cx(ngs,lc:lhab) 
+!     real ::  xv(ngs,lc:lhab) 
+!     real ::  vtxbar(ngs,lc:lhab,3) 
+!     real ::  xmas(ngs,lc:lhab) 
+!     real ::  xdn(ngs,lc:lhab) 
+!     real ::  xdia(ngs,lc:lhab,3) 
+!     real ::  vx(ngs,li:lhab) 
+!     real ::  alpha(ngs,lc:lhab) 
+!     real ::  zx(ngs,lr:lhab) 
+!     logical :: hasmass(nx,lc+1:lhab)
+!
+!     integer igs(ngs),kgs(ngs)
+!     
+!     real rho0(ngs),temcg(ngs)
+!
+!     real temg(ngs)
+!     
+!     real rhovt(ngs)
+!     
+!     real cwnc(ngs),cinc(ngs)
+!     real fadvisc(ngs),cwdia(ngs),cipmas(ngs)
+!     
+!     real cimasn,cimasx,cnina(ngs),cimas(ngs)
+!     
+!     real cnostmp(ngs)
+
+      real, allocatable ::  qx(:,:)
+      real, allocatable ::  qxw(:,:)
+      real, allocatable ::  cx(:,:)
+      real, allocatable ::  xv(:,:)
+      real, allocatable ::  vtxbar(:,:,:)
+      real, allocatable ::  xmas(:,:)
+      real, allocatable ::  xdn(:,:)
+      real, allocatable ::  xdia(:,:,:)
+      real, allocatable ::  vx(:,:)
+      real, allocatable ::  alpha(:,:)
+      real, allocatable ::  zx(:,:)
+      logical, allocatable :: hasmass(:,:)
+
+      integer, allocatable :: igs(:),kgs(:)
+      
+      real, allocatable :: rho0(:),temcg(:)
+
+      real, allocatable :: temg(:)
+      
+      real, allocatable :: rhovt(:)
+      
+      real, allocatable :: cwnc(:),cinc(:)
+      real, allocatable :: fadvisc(:),cwdia(:),cipmas(:)
+      
+      real, allocatable :: cnina(:),cimas(:)
+      
+      real, allocatable :: cnostmp(:)
+      
+      real :: cimasn,cimasx
       
 
 !-----------------------------------------------------------------------------
@@ -4336,7 +4370,30 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, &
 ! ###################################################################
 
 
-
+      allocate( db1(nx,nz+1),dtz1(nz+1,nx,0:1),dz2dinv(nz+1,nx),db1inv(nx,nz+1),rhovtzx(nz,nx) )
+      allocate( xfall0(nx,ny), xvt(nz+1,nx,3,lc:lhab), tmpn(-nor+ng1:nx+nor,-nor+ng1:ny+nor,-norz+ng1:nz+norz) )
+      allocate( tmpn2(-nor+ng1:nx+nor,-nor+ng1:ny+nor,-norz+ng1:nz+norz), z(-nor+ng1:nx+nor,-norz+ng1:nz+norz,lr:lhab))
+
+      ngs = nz+3
+      
+      allocate( qx(ngs,lv:lhab),  &
+                qxw(ngs,ls:lhab),  &
+                cx(ngs,lc:lhab),  &
+                xv(ngs,lc:lhab),  &
+                vtxbar(ngs,lc:lhab,3),  &
+                xmas(ngs,lc:lhab),  &
+                xdn(ngs,lc:lhab),  &
+                xdia(ngs,lc:lhab,3),  &
+                vx(ngs,li:lhab),  &
+                alpha(ngs,lc:lhab),  &
+                zx(ngs,lr:lhab),     &
+                hasmass(nx,lc+1:lhab), &
+                igs(ngs),kgs(ngs), &
+                rho0(ngs),temcg(ngs),temg(ngs), rhovt(ngs), &
+                cwnc(ngs),cinc(ngs), &
+                fadvisc(ngs),cwdia(ngs),cipmas(ngs), &
+                cnina(ngs),cimas(ngs), &
+                cnostmp(ngs) )
 
       kzb = 1
       kze = nz
@@ -4656,8 +4713,29 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, &
       ENDDO ! ix
 
 
+      deallocate( db1,dtz1,dz2dinv,db1inv,rhovtzx )
+      deallocate( xfall0, xvt, tmpn )
+      deallocate( tmpn2, z)
+
+      deallocate( qx,  &
+                qxw,  &
+                cx,  &
+                xv,  &
+                vtxbar,  &
+                xmas,  &
+                xdn,  &
+                xdia,  &
+                vx,  &
+                alpha,  &
+                zx,     &
+                hasmass, &
+                igs,kgs, &
+                rho0,temcg,temg, rhovt, &
+                cwnc,cinc, &
+                fadvisc,cwdia,cipmas, &
+                cnina,cimas, &
+                cnostmp )
 
-      
       RETURN
       END SUBROUTINE SEDIMENT1D