Skip to content

Commit

Permalink
Squashed commit of the following:
Browse files Browse the repository at this point in the history 
commit cdbde8389542f73be8f511b8c69d674639f82b0d
Author: Zhengui Wang <[email protected]>
Date:   Mon Jul 15 17:38:16 2024 -0400

    add more diagostic variables for ICM SAV module

commit 524d233b27f8e474e82752697a1cf62bf6513ddd
Author: Zhengui Wang <[email protected]>
Date:   Mon Jul 15 17:11:42 2024 -0400

    add more diagostic variables for ICM benthic algae

commit 142fc9bcdac76ef2780b27dfc2c240eb0879e021
Author: Zhengui Wang <[email protected]>
Date:   Mon Jul 15 15:44:57 2024 -0400

    more edit on ICM debug variables

commit 7bddfcf3ba7e03503f2df5f356bedebe514eadaa
Author: Zhengui Wang <[email protected]>
Date:   Mon Jul 15 15:30:56 2024 -0400

    rename the debug variables of ICM clam model

commit dd3ff74b6a892ffc930506c0b84f98e44e6ad848
Author: Zhengui Wang <[email protected]>
Date:   Thu Jul 11 18:23:56 2024 -0400

    1). work on CLAM diagostic variables; 2). fix bugs in the clam model

commit 877c02c469b3904ce0214070406e9f5311498ac8
Author: Zhengui Wang <[email protected]>
Date:   Fri Jul 5 12:31:11 2024 -0400

    add debug varaibles for ICM clam module
  • Loading branch information
@wzhengui
wzhengui committed Jul 15, 2024
1 parent 3074779 commit 83f51e6
Show file tree
Hide file tree
Showing 3 changed files with 113 additions and 30 deletions.
54 changes: 36 additions & 18 deletions src/ICM/icm.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -523,18 +523,19 @@ subroutine ecosystem(it)
!----------------------------------------------------------------------------------
if(iof_icm_dbg/=0) then
!Core
dbTN=0; dbTP=0 !TN, TP
db_TN=0; db_TP=0 !TN, TP
do k=kb+1,nvrt
dzb=(zid(k)-zid(k-1))
dbTN(id)=dzb*(RPON(k)+LPON(k)+DON(k)+NH4(k)+NO3(k)+sum(n2c*PBS(1:3,k)))
dbTP(id)=dzb*(RPOP(k)+LPOP(k)+DOP(k)+PO4(k)+sum(p2c*PBS(1:3,k)))
db_TN(id)=dzb*(RPON(k)+LPON(k)+DON(k)+NH4(k)+NO3(k)+sum(n2c*PBS(1:3,k)))
db_TP(id)=dzb*(RPOP(k)+LPOP(k)+DOP(k)+PO4(k)+sum(p2c*PBS(1:3,k)))
if(iCBP==1) then
dbTN(id)=dbTN(id)+dzb*SRPON(k)
dbTP(id)=dbTP(id)+dzb*SRPOP(k)
db_TN(id)=db_TN(id)+dzb*SRPON(k)
db_TP(id)=db_TP(id)+dzb*SRPOP(k)
endif
enddo
do k=kb+1,nvrt
dbCHLA(k,id)=max(sum(PBS(1:3,k)/c2chl(1:3)),0.d0) !CHLA
db_CHLA(k,id)=max(sum(PBS(1:3,k)/c2chl(1:3)),0.d0) !CHLA
db_Ke(k,id)=rKe(k)
enddo
endif !iof_icm_dbg/=0

Expand Down Expand Up @@ -859,7 +860,7 @@ subroutine zoo_calc(kb,PR)
end subroutine zoo_calc

subroutine sav_calc(id,kb,zid,shtz,tdep,sLight)
use schism_glbl, only : rkind,errmsg,idry_e,nvrt
use schism_glbl, only : rkind,errmsg,idry_e,nvrt,iof_icm_dbg
use icm_mod
implicit none
integer,intent(in) :: id,kb
Expand Down Expand Up @@ -911,6 +912,12 @@ subroutine sav_calc(id,kb,zid,shtz,tdep,sLight)
stem0=stem0-MT0(2)*dtw
sroot=sroot+(sFCP(3)*(1.0-sFAM)*sum(GP((kb+1):nvrt)*dz((kb+1):nvrt))-MT0(3))*dtw
endif

!debug
if(iof_icm_dbg/=0) then
db_sGP(k,id)=GP(k); db_sMT1(k,id)=MT1(k); db_sMT2(k,id)=MT2(k); db_sMT01(id)=MT0(1); db_sMT02(id)=MT0(2)
db_sMT03(id)=MT0(3); db_sfT=fT; db_sfI=fR; db_sfN=fN; db_sfP=fP
endif
enddo
sleaf=max(leaf0+sum(zleaf*dz),1.d-5)
sstem=max(stem0+sum(zstem*dz),1.d-5)
Expand Down Expand Up @@ -966,7 +973,7 @@ subroutine ba_calc(id,kb,wdz)
!----------------------------------------------------------------------------
!Benthic algae computation
!----------------------------------------------------------------------------
use schism_glbl,only : rkind
use schism_glbl,only : rkind,iof_icm_dbg
use icm_mod
implicit none
integer,intent(in) :: id,kb
Expand Down Expand Up @@ -1024,6 +1031,10 @@ subroutine ba_calc(id,kb,wdz)
JNH4(id)=JNH4(id)-gn2c*GP*fPN*(1.0-fWN)
JNO3(id)=JNO3(id)-gn2c*GP*(1.0-fPN)*(1.0-fWN)
JPO4(id)=JPO4(id)-gp2c*GP*(1.0-fWP)

if(iof_icm_dbg/=0) then
db_gfT(id)=fT; db_gfI(id)=fR; db_gfN(id)=fN; db_gfP(id)=fP
endif
endif !iBA==1

end subroutine ba_calc
Expand All @@ -1032,7 +1043,7 @@ subroutine clam_calc(id,kb,wdz)
!----------------------------------------------------------------------------
!clam model computation
!----------------------------------------------------------------------------
use schism_glbl,only : rkind,nvrt
use schism_glbl,only : rkind,nvrt,iof_icm_dbg
use schism_msgp, only : myrank,parallel_abort
use icm_mod
implicit none
Expand Down Expand Up @@ -1088,18 +1099,25 @@ subroutine clam_calc(id,kb,wdz)
if(idoy>=cDoyp(i,1).and.idoy<=cDoyp(i,2)) PR(i) =cPRR(i)*CLAM(i,id) !predation (g[C].m-2.day-1)
if(idoy>=cDoyh(i,1).and.idoy<=cDoyh(i,2)) HST(i)=cHSR(i)*CLAM(i,id) !harvest (g[C].m-2.day-1)
CLAM(i,id)=CLAM(i,id)+(GP(i)-MT(i)-RT(i)-PR(i)-HST(i))*dtw !update clam biomass

!debug
if(iof_icm_dbg/=0) then
db_cfT(i,id)=fT(i); db_cfS(i,id)=fS(i); db_cfDO(i,id)=fDO(i); db_cfTSS(i,id)=fTSS(i); db_cFr(i,id)=Fr(i)
db_cIF(i,id)=cIF(i); db_cTFC(i,id)=TFC(i); db_cATFC(i,id)=ATFC(i); db_cfN(i,id)=fN(i); db_cGP(i,id)=GP(i)
db_cMT(i,id)=MT(i); db_cRT(i,id)=RT(i); db_cPR(i,id)=PR(i); db_cHST(i,id) =HST(i)
endif
enddo !i=1,nclam

!interaction with water column variables; change rate of conc. (g.m-3.day-1)
cdwqc(iPB1, kb+1)=sum(PC(1)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iPB2, kb+1)=sum(PC(2)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iPB3, kb+1)=sum(PC(3)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iLPOC,kb+1)=sum(PC(4)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iRPOC,kb+1)=sum(PC(5)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iLPON,kb+1)=sum(PN(4)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iRPON,kb+1)=sum(PN(5)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iLPOP,kb+1)=sum(PP(4)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iRPOP,kb+1)=sum(PP(5)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iPB1, kb+1)=-sum(PC(1)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iPB2, kb+1)=-sum(PC(2)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iPB3, kb+1)=-sum(PC(3)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iLPOC,kb+1)=-sum(PC(4)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iRPOC,kb+1)=-sum(PC(5)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iLPON,kb+1)=-sum(PN(4)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iRPON,kb+1)=-sum(PN(5)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iLPOP,kb+1)=-sum(PP(4)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iRPOP,kb+1)=-sum(PP(5)*Fr*CLAM(1:nclam,id))/wdz
cdwqc(iNH4, kb+1)=sum((ATFN-cn2c*GP)+cn2c*MT)/wdz
cdwqc(iPO4, kb+1)=sum((ATFP-cp2c*GP)+cp2c*MT)/wdz
cdwqc(iDOX, kb+1)=o2c*sum((ATFC-GP)+MT)/wdz
Expand Down
76 changes: 67 additions & 9 deletions src/ICM/icm_init.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,7 +180,7 @@ subroutine read_icm_param(imode)
!p%itype is consistent with the output type in SCHISM hydro
!-----------------------------------------------------------------------------------------
!allocate variables
allocate(wqout(100),wqhot(50),stat=istat)
allocate(wqout(200),wqhot(50),stat=istat)
if(istat/=0) call parallel_abort('failed in alloc. wqout')

nouts=0; nout=0; iout=0; nhot_icm=0 !init
Expand All @@ -199,11 +199,13 @@ subroutine read_icm_param(imode)
nb=17; nouts(1)=nb; iout(1,1)=1; iout(2,1)=nb; nout=nout+nb

!debug variable
p=>wqout(nout+1); p%name='ICM_TN'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); dbTN=>p%p1; p%itype=4
p=>wqout(nout+2); p%name='ICM_TP'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); dbTP=>p%p1; p%itype=4
p=>wqout(nout+3); p%name='ICM_CHLA'; allocate(p%data(1,nvrt,nea)); p%p2=>p%data(1,:,:); dbCHLA=>p%p2; p%itype=6
!nb=3; nouts(1)=nouts(1)+nb; iout(2,1)=iout(2,1)+nb; nout=nout+nb
nb=3; nouts(1)=nouts(1)+nb; nout=nout+nb
if(iof_icm_dbg/=0) then
p=>wqout(nout+1); p%name='ICM_TN'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_TN=>p%p1; p%itype=4
p=>wqout(nout+2); p%name='ICM_TP'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_TP=>p%p1; p%itype=4
p=>wqout(nout+3); p%name='ICM_CHLA'; allocate(p%data(1,nvrt,nea)); p%p2=>p%data(1,:,:); db_CHLA=>p%p2; p%itype=6
p=>wqout(nout+4); p%name='ICM_Ke'; allocate(p%data(1,nvrt,nea)); p%p2=>p%data(1,:,:); db_Ke=>p%p2; p%itype=6
nb=4; nouts(1)=nouts(1)+nb; nout=nout+nb
endif
endif

!Silica module:2
Expand Down Expand Up @@ -266,6 +268,21 @@ subroutine read_icm_param(imode)
p=>wqout(nout+4); p%name='ICM_sht'; p%p1=>sht; p%itype=4
nb=4; nouts(6)=nb; iout(1,6)=nout+1; iout(2,6)=nout+nb; nout=nout+nb

!debug variable
if(iof_icm_dbg/=0) then
p=>wqout(nout+1); p%name='ICM_sGP'; allocate(p%data(1,nvrt,nea)); p%p2=>p%data(1,:,:); db_sGP=>p%p2; p%itype=6
p=>wqout(nout+2); p%name='ICM_sMT1'; allocate(p%data(1,nvrt,nea)); p%p2=>p%data(1,:,:); db_sMT1=>p%p2; p%itype=6
p=>wqout(nout+3); p%name='ICM_sMT2'; allocate(p%data(1,nvrt,nea)); p%p2=>p%data(1,:,:); db_sMT2=>p%p2; p%itype=6
p=>wqout(nout+4); p%name='ICM_sMT01'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_sMT01=>p%p1; p%itype=4
p=>wqout(nout+5); p%name='ICM_sMT02'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_sMT02=>p%p1; p%itype=4
p=>wqout(nout+6); p%name='ICM_sMT03'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_sMT03=>p%p1; p%itype=4
p=>wqout(nout+7); p%name='ICM_sfT'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_sfT=>p%p1; p%itype=4
p=>wqout(nout+8); p%name='ICM_sfI'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_sfI=>p%p1; p%itype=4
p=>wqout(nout+9); p%name='ICM_sfN'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_sfN=>p%p1; p%itype=4
p=>wqout(nout+10); p%name='ICM_sfP'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_sfP=>p%p1; p%itype=4
nb=10; nouts(1)=nouts(1)+nb; nout=nout+nb
endif

!hotstart variable
p=>wqhot(nhot+1); p%name='sav'; p%p2=>sav
nhot=nhot+1
Expand Down Expand Up @@ -342,11 +359,20 @@ subroutine read_icm_param(imode)

!Benthic Algea: 9
if(iBA==1) then
p=>wqout(nout+1); p%name='ICM_gBA'; p%p1=>gBA; p%itype=4
p=>wqout(nout+1); p%name='ICM_gBA'; p%p1=>gBA; p%itype=4
p=>wqout(nout+2); p%name='ICM_gGP'; p%p1=>gGP; p%itype=4
p=>wqout(nout+3); p%name='ICM_gMT'; p%p1=>gMT; p%itype=4
p=>wqout(nout+4); p%name='ICM_gPR'; p%p1=>gPR; p%itype=4
nb=4; nouts(9)=nb; iout(1,9)=nout+1; iout(2,9)=nout+nb; nout=nout+nb

!debug variable
if(iof_icm_dbg/=0) then
p=>wqout(nout+1); p%name='ICM_gfT'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_gfT=>p%p1; p%itype=4
p=>wqout(nout+2); p%name='ICM_gfI'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_gfI=>p%p1; p%itype=4
p=>wqout(nout+3); p%name='ICM_gfN'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_gfN=>p%p1; p%itype=4
p=>wqout(nout+4); p%name='ICM_gfP'; allocate(p%data(1,1,nea)); p%p1=>p%data(1,1,:); db_gfP=>p%p1; p%itype=4
nb=4; nouts(1)=nouts(1)+nb; nout=nout+nb
endif

!hotstart
p=>wqhot(nhot+1); p%name='gBA'; p%p1=>gBA
Expand All @@ -361,6 +387,28 @@ subroutine read_icm_param(imode)
enddo
nb=nclam; nouts(10)=nb; iout(1,10)=nout+1; iout(2,10)=nout+nb; nout=nout+nb

!debug variable
if(iof_icm_dbg/=0) then
do i=1,nclam
write(stmp,"(I3)") i
p=>wqout(nout+(i-1)*14+1); p%name='ICM_cfT'//trim(adjustl(stmp)); p%p1=>db_cfT(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+2); p%name='ICM_cfS'//trim(adjustl(stmp)); p%p1=>db_cfS(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+3); p%name='ICM_cfDO'//trim(adjustl(stmp)); p%p1=>db_cfDO(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+4); p%name='ICM_cfTSS'//trim(adjustl(stmp)); p%p1=>db_cfTSS(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+5); p%name='ICM_cFr'//trim(adjustl(stmp)); p%p1=>db_cFr(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+6); p%name='ICM_cIF'//trim(adjustl(stmp)); p%p1=>db_cIF(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+7); p%name='ICM_cTFC'//trim(adjustl(stmp)); p%p1=>db_cTFC(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+8); p%name='ICM_cATFC'//trim(adjustl(stmp)); p%p1=>db_cATFC(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+9); p%name='ICM_cfN'//trim(adjustl(stmp)); p%p1=>db_cfN(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+10); p%name='ICM_cGP'//trim(adjustl(stmp)); p%p1=>db_cGP(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+11); p%name='ICM_cMT'//trim(adjustl(stmp)); p%p1=>db_cMT(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+12); p%name='ICM_cRT'//trim(adjustl(stmp)); p%p1=>db_cRT(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+13); p%name='ICM_cPR'//trim(adjustl(stmp)); p%p1=>db_cPR(i,:); p%itype=4
p=>wqout(nout+(i-1)*14+14); p%name='ICM_cHST'//trim(adjustl(stmp)); p%p1=>db_cHST(i,:); p%itype=4
enddo
nb=14*nclam; nouts(1)=nouts(1)+nb; nout=nout+nb
endif

!hotstart
p=>wqhot(nhot+1); p%name='clam'; p%p2=>CLAM
nhot=nhot+1
Expand Down Expand Up @@ -393,7 +441,7 @@ subroutine read_icm_param(imode)
elseif(associated(p%p3)) then
p%ndim=3; p%dims(1:3)=shape(p%p3)
endif
if(p%dims(3)/=nea) call parallel_abort('ICM hotstart variable dim/=nea')
if(p%dims(3)/=nea) call parallel_abort('ICM hotstart variable dim/=nea: '//trim(adjustl(p%name)))
enddo

!------------------------------------------------------------------------------------
Expand Down Expand Up @@ -595,7 +643,8 @@ subroutine icm_vars_init
!allocate ICM variables and initialize; read ICM spatial parameters
!--------------------------------------------------------------------------------
use schism_glbl, only : rkind,nea,npa,nvrt,irange_tr,ntrs,in_dir,len_in_dir,np_global, &
& ne_global,ielg,iplg,i34,elnode,flag_ic,tr_nd,tr_el,indel,np,nne
& ne_global,ielg,iplg,i34,elnode,flag_ic,tr_nd,tr_el,indel,np, &
& nne,iof_icm_dbg
use schism_msgp, only : exchange_p3d_tr,parallel_abort,myrank,comm,itype,rtype
use netcdf
use icm_mod
Expand Down Expand Up @@ -671,6 +720,15 @@ subroutine icm_vars_init
if(iClam==1) then
allocate(CLAM(nclam,nea),cFPOC(nea,2),cFPON(nea,2),cFPOP(nea,2), stat=istat)
if(istat/=0) call parallel_abort('failed in alloc. CLAM')
if(iof_icm_dbg/=0) then
allocate(db_cfT(nclam,nea),db_cfS(nclam,nea),db_cfDO(nclam,nea),db_cfTSS(nclam,nea), &
& db_cFr(nclam,nea),db_cIF(nclam,nea),db_cTFC(nclam,nea),db_cATFC(nclam,nea),&
& db_cfN(nclam,nea),db_cGP(nclam,nea),db_cMT(nclam,nea), db_cRT(nclam,nea), &
& db_cPR(nclam,nea),db_cHST(nclam,nea), stat=istat)
if(istat/=0) call parallel_abort('failed in alloc. db_cfT')
db_cfT=0; db_cfS=0; db_cfDO=0; db_cfTSS=0; db_cFr=0; db_cIF=0; db_cTFC=0; db_cATFC=0
db_cfN=0; db_cGP=0; db_cMT=0; db_cRT=0; db_cPR=0; db_cHST=0
endif
endif

!-------------------------------------------------------------------------------
Expand Down
13 changes: 10 additions & 3 deletions src/ICM/icm_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,8 +53,8 @@ module icm_mod
real(rkind),save,pointer,dimension(:) :: zdDOX

!debug variables
real(rkind),save,pointer,dimension(:,:) :: dbCHLA
real(rkind),save,pointer,dimension(:) :: dbTN,dbTP
real(rkind),save,pointer,dimension(:,:) :: db_CHLA,db_Ke
real(rkind),save,pointer,dimension(:) :: db_TN,db_TP

!-------------------------------------------------------------------------------
!ICM parameters and variables
Expand Down Expand Up @@ -110,6 +110,8 @@ module icm_mod

integer,save,allocatable :: spatch(:) !sav region
real(rkind),save,target,allocatable :: sht(:),sav(:,:),sFPOC(:,:),sFPON(:,:),sFPOP(:,:),sbNH4(:),sbPO4(:),sSOD(:)
real(rkind),save,pointer,dimension(:,:) :: db_sGP,db_sMT1,db_sMT2
real(rkind),save,pointer,dimension(:) :: db_sMT01,db_sMT02,db_sMT03,db_sfT,db_sfI,db_sfN,db_sfP

!-------------------------------------------------------------------------------
!marsh parameters and variables
Expand Down Expand Up @@ -160,7 +162,7 @@ module icm_mod

integer,save,allocatable,dimension(:) :: gpatch
real(rkind),save,target,allocatable,dimension(:) :: gBA,gGP,gMT,gPR

real(rkind),save,pointer,dimension(:) :: db_gfT,db_gfI,db_gfN,db_gfP
!-------------------------------------------------------------------------------
!Clam model (CLAM) parameters and variables
!-------------------------------------------------------------------------------
Expand All @@ -171,6 +173,11 @@ module icm_mod

integer,save,allocatable,dimension(:) :: cpatch
real(rkind),save,target,allocatable,dimension(:,:) :: CLAM,cFPOC,cFPON,cFPOP

!debug
real(rkind),save,pointer,dimension(:,:) :: db_cfT,db_cfS,db_cfDO,db_cfTSS,db_cFr,db_cIF,db_cTFC,db_cATFC,db_cfN, &
& db_cGP,db_cMT,db_cRT,db_cPR,db_cHST

!-------------------------------------------------------------------------------
!benthic erosion (ERO) parameters and variables
!-------------------------------------------------------------------------------
Expand Down

0 comments on commit 83f51e6

Please sign in to comment.