From f416bb47d2057612d6212883b34626718d4220f9 Mon Sep 17 00:00:00 2001
From: Benoit Urruty <benoit.urruty@univ-grenoble-alpes.fr>
Date: Tue, 9 Feb 2021 19:57:11 +0100
Subject: [PATCH] algo
---
src/PICO.F90 | 1013 ++++++++++++++++++++--------------------
src/boxmodel_solver.md | 296 +++++++++++-
2 files changed, 801 insertions(+), 508 deletions(-)
diff --git a/src/PICO.F90 b/src/PICO.F90
index 0f7f9af..c1c0ef2 100644
--- a/src/PICO.F90
+++ b/src/PICO.F90
@@ -1,555 +1,556 @@
SUBROUTINE PICO( Model,Solver,dt,Transient )
-!------------------------------------------------------------------------------
-!USE CoordinateSystems
-!USE MeshUtils
-USE Netcdf
-USE DefUtils
-
-IMPLICIT NONE
-!------------------------------------------------------------------------------
-TYPE(Model_t) :: Model
-TYPE(Solver_t), TARGET :: Solver
-LOGICAL :: Transient
-REAL(KIND=dp) :: dt
-
-!------------------------------------------------------------------------------
-! Local variables
-!------------------------------------------------------------------------------
-TYPE(Mesh_t),POINTER :: Mesh
-TYPE(ValueList_t), POINTER :: Params
-TYPE(Solver_t),POINTER :: PSolver
-TYPE(Variable_t),POINTER :: MeltVar=>NULL(), GMVar=>NULL(), DepthVar=>NULL(), SVar=>NULL(), TVar=>NULL(), BoxVar=>NULL()
-TYPE(Variable_t),POINTER :: isfslopeVar=>NULL(), distGLVar=>NULL(), distIFVar=>NULL()
-TYPE(Variable_t),POINTER :: TimeVar=>NULL()
-TYPE(Nodes_t) :: ElementNodes
-TYPE(GaussIntegrationPoints_t) :: IntegStuff
-TYPE(Element_t),POINTER :: Element
-
-
-REAL(kind=dp),ALLOCATABLE :: VisitedNode(:),db(:),Basis(:),dBasisdx(:,:) ,Depth(:)
-REAL(kind=dp) :: u,v,w,SqrtElementMetric,s
-
-INTEGER , POINTER :: MeltPerm(:), GMPerm(:), DepthPerm(:),NodeIndexes(:), SPerm(:), TPerm(:), BPerm(:), loc(:), isfslopePerm(:), distGLPerm(:), distIFPerm(:)
-INTEGER , DIMENSION(:), ALLOCATABLE :: boxes(:)
-REAL(KIND=dp) , POINTER :: Melt(:),GM(:),isfslope(:), Boxnumber(:), DATAPointer(:,:), distGL(:), distIF(:),DepthVal(:)
-
-LOGICAL :: stat, Found,UnFoundFatal=.TRUE.
-
-CHARACTER(len=MAX_NAME_LEN) :: variabletype, VariableName
-CHARACTER(len = 200) :: meltValue
-CHARACTER(LEN=MAX_NAME_LEN) :: SolverName='MELT_MISMIP', FName, DataF
-
-INTEGER :: node, e, t, n, i, j, knd, kk, nD, ii, b, ierr, NetcdfStatus,varid,ncid
-INTEGER :: status(MPI_STATUS_SIZE)
-INTEGER :: maxbastmp
-
-REAL(KIND=dp) :: localInteg, Integ, Integ_Reduced, zzz, tmp1, xbox, Tstar, &
-& localunity, Area_Reduced, g1, g2, rr,nn, &
-& sn, distmax, max_Reduced, TT, SS, qqq_Reduced
-REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: basinmax,basin_Reduced , S0, T0,qqq
-REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: Mbox, Tbox, Sbox, Zbox, Abox
+ !------------------------------------------------------------------------------
+ !USE CoordinateSystems
+ !USE MeshUtils
+ USE Netcdf
+ USE DefUtils
+
+ IMPLICIT NONE
+ !------------------------------------------------------------------------------
+ TYPE(Model_t) :: Model
+ TYPE(Solver_t), TARGET :: Solver
+ LOGICAL :: Transient
+ REAL(KIND=dp) :: dt
+
+ !------------------------------------------------------------------------------
+ ! Local variables
+ !------------------------------------------------------------------------------
+ TYPE(Mesh_t),POINTER :: Mesh
+ TYPE(ValueList_t), POINTER :: Params
+ TYPE(Solver_t),POINTER :: PSolver
+ TYPE(Variable_t),POINTER :: MeltVar=>NULL(), GMVar=>NULL(), DepthVar=>NULL(), SVar=>NULL(), TVar=>NULL(), BoxVar=>NULL()
+ TYPE(Variable_t),POINTER :: isfslopeVar=>NULL(), distGLVar=>NULL(), distIFVar=>NULL()
+ TYPE(Variable_t),POINTER :: TimeVar=>NULL()
+ TYPE(Nodes_t) :: ElementNodes
+ TYPE(GaussIntegrationPoints_t) :: IntegStuff
+ TYPE(Element_t),POINTER :: Element
+
+
+ REAL(kind=dp),ALLOCATABLE :: VisitedNode(:),db(:),Basis(:),dBasisdx(:,:) ,Depth(:)
+ REAL(kind=dp) :: u,v,w,SqrtElementMetric,s
+
+ INTEGER , POINTER :: MeltPerm(:), GMPerm(:), DepthPerm(:),NodeIndexes(:), SPerm(:), TPerm(:), BPerm(:), loc(:), isfslopePerm(:), distGLPerm(:), distIFPerm(:)
+ INTEGER , DIMENSION(:), ALLOCATABLE :: boxes(:)
+ REAL(KIND=dp) , POINTER :: Melt(:),GM(:),isfslope(:), Boxnumber(:), DATAPointer(:,:), distGL(:), distIF(:),DepthVal(:)
+
+ LOGICAL :: stat, Found,UnFoundFatal=.TRUE.
+
+ CHARACTER(len=MAX_NAME_LEN) :: variabletype, VariableName
+ CHARACTER(len = 200) :: meltValue
+ CHARACTER(LEN=MAX_NAME_LEN) :: SolverName='MELT_MISMIP', FName, DataF
+
+ INTEGER :: node, e, t, n, i, j, knd, kk, nD, ii, b, ierr, NetcdfStatus,varid,ncid
+ INTEGER :: status(MPI_STATUS_SIZE)
+ INTEGER :: maxbastmp
+
+ REAL(KIND=dp) :: localInteg, Integ, Integ_Reduced, zzz, tmp1, xbox, Tstar, &
+ & localunity, Area_Reduced, g1, g2, rr,nn, &
+ & sn, distmax, max_Reduced, TT, SS, qqq_Reduced
+ REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: basinmax,basin_Reduced , S0, T0,qqq
+ REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: Mbox, Tbox, Sbox, Zbox, Abox
!!!! SAVE
-TYPE(Variable_t),POINTER, SAVE :: BasinVar=>NULL()
-LOGICAL,SAVE :: Initialized = .FALSE., ExtrudedMesh=.False.
-LOGICAL, SAVE :: Firsttime=.TRUE., llGL, Parallel
-INTEGER, SAVE :: boxmax, Nmax, MaxBas
-INTEGER :: tmeanid, nlen
-INTEGER, POINTER, SAVE :: BasinPerm(:)
-REAL(KIND=dp), POINTER, SAVE :: Basin(:)
-REAL(KIND=dp), DIMENSION(:), ALLOCATABLE, SAVE :: S_mean, T_mean
-REAL(KIND=dp), SAVE :: sealevel, lbd1, lbd2, lbd3, meltfac, K, gT, rhostar, CC,beta, alpha, mskcrit
+ TYPE(Variable_t),POINTER, SAVE :: BasinVar=>NULL()
+ LOGICAL,SAVE :: Initialized = .FALSE., ExtrudedMesh=.False.
+ LOGICAL, SAVE :: Firsttime=.TRUE., llGL, Parallel
+ INTEGER, SAVE :: boxmax, Nmax, MaxBas
+ INTEGER :: tmeanid, nlen
+ INTEGER, POINTER, SAVE :: BasinPerm(:)
+ REAL(KIND=dp), POINTER, SAVE :: Basin(:)
+ REAL(KIND=dp), DIMENSION(:), ALLOCATABLE, SAVE :: S_mean, T_mean
+ REAL(KIND=dp), SAVE :: sealevel, lbd1, lbd2, lbd3, meltfac, K, gT, rhostar, CC,beta, alpha, mskcrit
-!------------------------------------------------------------------------------
-! 1- Read constants and parameters of the simulation :
-!------------------------------------------------------------------------------
+ !------------------------------------------------------------------------------
+ ! 1- Read constants and parameters of the simulation :
+ !------------------------------------------------------------------------------
-!- Simulation parameters (idealized ocean forcing) :
+ !- Simulation parameters (idealized ocean forcing) :
-Params => GetSolverParams()
+ Params => GetSolverParams()
-Mesh => Model % Mesh
+ Mesh => Model % Mesh
-Nmax = Solver % Mesh % NumberOfNodes
-BasinVar => VariableGet( Model % Mesh % Variables, 'basins')
-IF (.NOT.ASSOCIATED(BasinVar)) &
- & CALL FATAL(SolverName,'basins not found')
-print *,'Vartype',BasinVar % TYPE
+ Nmax = Solver % Mesh % NumberOfNodes
+ BasinVar => VariableGet( Model % Mesh % Variables, 'basins')
+ IF (.NOT.ASSOCIATED(BasinVar)) &
+ & CALL FATAL(SolverName,'basins not found')
+ print *,'Vartype',BasinVar % TYPE
-!IF ( BasinVar % TYPE .NE. Variable_on_element) &
-! & CALL FATAL(SolverName,'basins not a vairable on element')
+ !IF ( BasinVar % TYPE .NE. Variable_on_element) &
+ ! & CALL FATAL(SolverName,'basins not a vairable on element')
-BasinPerm => BasinVar % Perm
-Basin => BasinVar % Values
+ BasinPerm => BasinVar % Perm
+ Basin => BasinVar % Values
-IF (Firsttime) THEN
- Firsttime=.False.
+ IF (Firsttime) THEN
+ Firsttime=.False.
- llGL=ListGetLogical( Model % Simulation, 'Grounding Line Melt', UnFoundFatal=UnFoundFatal )
-
- !- General :
+ llGL=ListGetLogical( Model % Simulation, 'Grounding Line Melt', UnFoundFatal=UnFoundFatal )
+
+ !- General :
+
+ sealevel = ListGetCReal( Model % Constants, 'Sea Level',UnFoundFatal=UnFoundFatal)
- sealevel = ListGetCReal( Model % Constants, 'Sea Level',UnFoundFatal=UnFoundFatal)
+ meltfac=ListGetCReal( Model % Constants, 'Melt factor',UnFoundFatal=UnFoundFatal)
- meltfac=ListGetCReal( Model % Constants, 'Melt factor',UnFoundFatal=UnFoundFatal)
+ lbd1 = ListGetCReal( Model % Constants, 'Liquidus slope',UnFoundFatal=UnFoundFatal )
- lbd1 = ListGetCReal( Model % Constants, 'Liquidus slope',UnFoundFatal=UnFoundFatal )
+ lbd2 = ListGetCReal( Model % Constants, 'Liquidus intercept',UnFoundFatal=UnFoundFatal)
- lbd2 = ListGetCReal( Model % Constants, 'Liquidus intercept',UnFoundFatal=UnFoundFatal)
+ lbd3 = ListGetCReal( Model % Constants, 'Liquidus pressure coeff',UnFoundFatal=UnFoundFatal)
- lbd3 = ListGetCReal( Model % Constants, 'Liquidus pressure coeff',UnFoundFatal=UnFoundFatal)
+ Parallel = (ParEnv %PEs > 1)
- Parallel = (ParEnv %PEs > 1)
+ maxbastmp=MAXVAL(Basin)
- maxbastmp=MAXVAL(Basin)
-
- IF (Parallel) THEN
- CALL MPI_ALLREDUCE(maxbastmp,MaxBas,1,MPI_INTEGER,MPI_MAX,ELMER_COMM_WORLD,ierr)
- ELSE
- MaxBas=maxbastmp
- END IF
+ IF (Parallel) THEN
+ CALL MPI_ALLREDUCE(maxbastmp,MaxBas,1,MPI_INTEGER,MPI_MAX,ELMER_COMM_WORLD,ierr)
+ ELSE
+ MaxBas=maxbastmp
+ END IF
- !ALLOCATE(S_mean(MaxBas), T_mean(MaxBas))
+ !ALLOCATE(S_mean(MaxBas), T_mean(MaxBas))
- DataF = ListGetString( Params, 'Data File', Found, UnFoundFatal )
+ DataF = ListGetString( Params, 'Data File', Found, UnFoundFatal )
- NetCDFstatus = NF90_OPEN(DataF,NF90_NOWRITE,ncid)
+ NetCDFstatus = NF90_OPEN(DataF,NF90_NOWRITE,ncid)
- NetCDFstatus = nf90_inq_dimid(ncid,'number_of_basins' , tmeanid)
- NetCDFstatus = nf90_inquire_dimension(ncid, tmeanid , len=nlen)
+ NetCDFstatus = nf90_inq_dimid(ncid,'number_of_basins' , tmeanid)
+ NetCDFstatus = nf90_inquire_dimension(ncid, tmeanid , len=nlen)
- ALLOCATE(S_mean(nlen), T_mean(nlen))
+ ALLOCATE(S_mean(nlen), T_mean(nlen))
- NetCDFstatus = nf90_inq_varid(ncid,'T_mean',varid)
- NetCDFstatus = nf90_get_var(ncid, varid,T_mean)
- IF ( NetCDFstatus /= NF90_NOERR ) THEN
- CALL Fatal(Trim(SolverName), &
- 'Unable to get netcdf variable T_mean')
- END IF
+ NetCDFstatus = nf90_inq_varid(ncid,'T_mean',varid)
+ NetCDFstatus = nf90_get_var(ncid, varid,T_mean)
+ IF ( NetCDFstatus /= NF90_NOERR ) THEN
+ CALL Fatal(Trim(SolverName), &
+ 'Unable to get netcdf variable T_mean')
+ END IF
- NetCDFstatus = nf90_inq_varid(ncid,'S_mean',varid)
- NetCDFstatus = nf90_get_var(ncid, varid,S_mean)
- IF ( NetCDFstatus /= NF90_NOERR ) THEN
- CALL Fatal(Trim(SolverName), &
- 'Unable to get netcdf variable S_mean')
- END IF
+ NetCDFstatus = nf90_inq_varid(ncid,'S_mean',varid)
+ NetCDFstatus = nf90_get_var(ncid, varid,S_mean)
+ IF ( NetCDFstatus /= NF90_NOERR ) THEN
+ CALL Fatal(Trim(SolverName), &
+ 'Unable to get netcdf variable S_mean')
+ END IF
- NetCDFstatus = nf90_inq_varid(ncid,'max box',varid)
- NetCDFstatus = nf90_get_var(ncid, varid,boxmax)
- IF ( NetCDFstatus /= NF90_NOERR ) THEN
- CALL Fatal(Trim(SolverName), &
- 'Unable to get netcdf max box')
- END IF
+ NetCDFstatus = nf90_inq_varid(ncid,'max box',varid)
+ NetCDFstatus = nf90_get_var(ncid, varid,boxmax)
+ IF ( NetCDFstatus /= NF90_NOERR ) THEN
+ CALL Fatal(Trim(SolverName), &
+ 'Unable to get netcdf max box')
+ END IF
- NetCDFstatus = nf90_inq_varid(ncid,'Circulation_Parameter',varid)
- NetCDFstatus = nf90_get_var(ncid, varid,CC)
- IF ( NetCDFstatus /= NF90_NOERR ) THEN
- CALL Fatal(Trim(SolverName), &
- 'Unable to get netcdf Circulation_Parameter')
- END IF
+ NetCDFstatus = nf90_inq_varid(ncid,'Circulation_Parameter',varid)
+ NetCDFstatus = nf90_get_var(ncid, varid,CC)
+ IF ( NetCDFstatus /= NF90_NOERR ) THEN
+ CALL Fatal(Trim(SolverName), &
+ 'Unable to get netcdf Circulation_Parameter')
+ END IF
- NetCDFstatus = nf90_inq_varid(ncid,'Effective Exchange Velocity',varid)
- NetCDFstatus = nf90_get_var(ncid, varid,gT)
- IF ( NetCDFstatus /= NF90_NOERR ) THEN
- CALL Fatal(Trim(SolverName), &
- 'Unable to get netcdf Effective Exchange Velocity')
- END IF
-
- NetCDFstatus = nf90_inq_varid(ncid,'Thermal Expansion coeff',varid)
- NetCDFstatus = nf90_get_var(ncid, varid,alpha)
- IF ( NetCDFstatus /= NF90_NOERR ) THEN
- CALL Fatal(Trim(SolverName), &
- 'Unable to get netcdf Thermal Expansion coeff')
- END IF
+ NetCDFstatus = nf90_inq_varid(ncid,'Effective Exchange Velocity',varid)
+ NetCDFstatus = nf90_get_var(ncid, varid,gT)
+ IF ( NetCDFstatus /= NF90_NOERR ) THEN
+ CALL Fatal(Trim(SolverName), &
+ 'Unable to get netcdf Effective Exchange Velocity')
+ END IF
- NetCDFstatus = nf90_inq_varid(ncid,'Salinity contraction coeff',varid)
- NetCDFstatus = nf90_get_var(ncid, varid,beta)
- IF ( NetCDFstatus /= NF90_NOERR ) THEN
- CALL Fatal(Trim(SolverName), &
- 'Unable to get netcdf Salinity contraction coeff')
- END IF
+ NetCDFstatus = nf90_inq_varid(ncid,'Thermal Expansion coeff',varid)
+ NetCDFstatus = nf90_get_var(ncid, varid,alpha)
+ IF ( NetCDFstatus /= NF90_NOERR ) THEN
+ CALL Fatal(Trim(SolverName), &
+ 'Unable to get netcdf Thermal Expansion coeff')
+ END IF
- NetCDFstatus = nf90_inq_varid(ncid,'EOS ref Density',varid)
- NetCDFstatus = nf90_get_var(ncid, varid,rhostar)
- IF ( NetCDFstatus /= NF90_NOERR ) THEN
- CALL Fatal(Trim(SolverName), &
- 'Unable to get netcdfEOS ref Density')
- END IF
- NetCDFstatus = nf90_close(ncid)
+ NetCDFstatus = nf90_inq_varid(ncid,'Salinity contraction coeff',varid)
+ NetCDFstatus = nf90_get_var(ncid, varid,beta)
+ IF ( NetCDFstatus /= NF90_NOERR ) THEN
+ CALL Fatal(Trim(SolverName), &
+ 'Unable to get netcdf Salinity contraction coeff')
+ END IF
- IF ( llGL ) THEN
- mskcrit = 0.5 ! Melt is at the Grounding Line and floating points
- else
- mskcrit = -0.5 ! No melt at the Grounding Line, only floating points
- ENDif
- CALL INFO(Trim(SolverName),'END FIRST TIME', Level =5)
-END IF
-CALL INFO(Trim(SolverName),'START', Level =5)
+ NetCDFstatus = nf90_inq_varid(ncid,'EOS ref Density',varid)
+ NetCDFstatus = nf90_get_var(ncid, varid,rhostar)
+ IF ( NetCDFstatus /= NF90_NOERR ) THEN
+ CALL Fatal(Trim(SolverName), &
+ 'Unable to get netcdf EOS ref Density')
+ END IF
+ NetCDFstatus = nf90_close(ncid)
-BoxVar => VariableGet( Model % Mesh % Variables, 'Boxes')
-IF (.NOT.ASSOCIATED(BoxVar)) &
- & CALL FATAL(SolverName,'boxes not found')
+ IF ( llGL ) THEN
+ mskcrit = 0.5 ! Melt is at the Grounding Line and floating points
+ else
+ mskcrit = -0.5 ! No melt at the Grounding Line, only floating points
+ ENDif
+ CALL INFO(Trim(SolverName),'END FIRST TIME', Level =5)
+ END IF
+ CALL INFO(Trim(SolverName),'START', Level =5)
-MeltVar => VariableGet( Model % Mesh % Variables, 'Melt')
-IF (.NOT.ASSOCIATED(MeltVar)) &
- & CALL FATAL(SolverName,'Melt not found')
-
-GMVar => VariableGet( Model % Mesh % Variables, 'GroundedMask')
-IF (.NOT.ASSOCIATED(GMVar)) &
- & CALL FATAL(SolverName,'GroundedMask not found')
-
-DepthVar => VariableGet( Model % Mesh % Variables, 'Zb')
-IF (.NOT.ASSOCIATED(DepthVar)) &
- & CALL FATAL(SolverName,'Zb not found')
-
-distGLVar => VariableGet( Model % Mesh % Variables, 'distGL')
-IF (.NOT.ASSOCIATED(distGLVar)) &
- & CALL FATAL(SolverName,'distGL not found')
-
-distIFVar => VariableGet( Model % Mesh % Variables, 'distIF')
-IF (.NOT.ASSOCIATED(distIFVar)) &
- & CALL FATAL(SolverName,'distIF not found')
-
-BPerm => BoxVar % Perm
-Boxnumber => BoxVar % Values
-
-MeltPerm => MeltVar % Perm
-Melt => MeltVar % Values
-
-GMPerm => GMVar % Perm
-GM => GMVar % Values
-
-DepthPerm => DepthVar % Perm
-DepthVal => DepthVar % Values
-Allocate (Depth(SIZE(DepthVal)))
-Depth = sealevel - DepthVal ! Depth < 0 under sea level
-
-ALLOCATE( Zbox(boxmax,MaxBas), Abox(boxmax,MaxBas), Tbox(boxmax,MaxBas), Sbox(boxmax,MaxBas), Mbox(boxmax,MaxBas),qqq(MaxBas), T0(MaxBas), S0(MaxBas))
-ALLOCATE(basin_Reduced(MaxBas),basinmax(MaxBas),boxes(MaxBas))
-ALLOCATE(VisitedNode(Nmax), Basis(Model % MaxElementNodes), dBasisdx(Model % MaxElementNodes,3))
-
-CALL INFO(Trim(SolverName),'LOAD Variables', Level =5)
-
-
-!test parameters, have to be change!!!!
-!dTemp = 0.0 ! warming (over time scale timsc)
-!timsc = 1.e9! time scale for changes dTempt (in years)
-
-T0 = T_mean
-S0 = S_mean
-
-distGLPerm => distGLVar % Perm
-distGL => distGLVar % Values
-
-distIFPerm => distIFVar % Perm
-distIF => distIFVar % Values
-
-VisitedNode=0.0_dp
-
-Abox(:,:)=0.0_dp
-Zbox(:,:)=0.0_dp
-basinmax(:)=0.0_dp
-distmax = 0.0_dp
-basin_Reduced(:)=0.0_dp
-boxes(:)=0.0_dp
+ BoxVar => VariableGet( Model % Mesh % Variables, 'Boxes')
+ IF (.NOT.ASSOCIATED(BoxVar)) &
+ & CALL FATAL(SolverName,'boxes not found')
+ MeltVar => VariableGet( Model % Mesh % Variables, 'Melt')
+ IF (.NOT.ASSOCIATED(MeltVar)) &
+ & CALL FATAL(SolverName,'Melt not found')
+
+ GMVar => VariableGet( Model % Mesh % Variables, 'GroundedMask')
+ IF (.NOT.ASSOCIATED(GMVar)) &
+ & CALL FATAL(SolverName,'GroundedMask not found')
+
+ DepthVar => VariableGet( Model % Mesh % Variables, 'Zb')
+ IF (.NOT.ASSOCIATED(DepthVar)) &
+ & CALL FATAL(SolverName,'Zb not found')
+
+ distGLVar => VariableGet( Model % Mesh % Variables, 'distGL')
+ IF (.NOT.ASSOCIATED(distGLVar)) &
+ & CALL FATAL(SolverName,'distGL not found')
+
+ distIFVar => VariableGet( Model % Mesh % Variables, 'distIF')
+ IF (.NOT.ASSOCIATED(distIFVar)) &
+ & CALL FATAL(SolverName,'distIF not found')
+
+ BPerm => BoxVar % Perm
+ Boxnumber => BoxVar % Values
+
+ MeltPerm => MeltVar % Perm
+ Melt => MeltVar % Values
+
+ GMPerm => GMVar % Perm
+ GM => GMVar % Values
+
+ DepthPerm => DepthVar % Perm
+ DepthVal => DepthVar % Values
+ Allocate (Depth(SIZE(DepthVal)))
+ Depth = sealevel - DepthVal ! Depth < 0 under sea level
+
+ ALLOCATE( Zbox(boxmax,MaxBas), Abox(boxmax,MaxBas), Tbox(boxmax,MaxBas), Sbox(boxmax,MaxBas), Mbox(boxmax,MaxBas),qqq(MaxBas), T0(MaxBas), S0(MaxBas))
+ ALLOCATE(basin_Reduced(MaxBas),basinmax(MaxBas),boxes(MaxBas))
+ ALLOCATE(VisitedNode(Nmax), Basis(Model % MaxElementNodes), dBasisdx(Model % MaxElementNodes,3))
+
+ CALL INFO(Trim(SolverName),'LOAD Variables', Level =5)
+
+
+ !test parameters, have to be change!!!!
+ !dTemp = 0.0 ! warming (over time scale timsc)
+ !timsc = 1.e9! time scale for changes dTempt (in years)
+
+ T0 = T_mean
+ S0 = S_mean
+
+ distGLPerm => distGLVar % Perm
+ distGL => distGLVar % Values
+
+ distIFPerm => distIFVar % Perm
+ distIF => distIFVar % Values
+
+ VisitedNode=0.0_dp
+
+ Abox(:,:)=0.0_dp
+ Zbox(:,:)=0.0_dp
+ basinmax(:)=0.0_dp
+ distmax = 0.0_dp
+ basin_Reduced(:)=0.0_dp
+ boxes(:)=0.0_dp
+ melt(:)=0.0_dp
CALL INFO(Trim(SolverName),'START BOXES', Level =5)
- ! first loop on element to determine the number of boxes per basins
-DO e=1,Solver % NumberOfActiveElements
- Element => GetActiveElement(e)
- CALL GetElementNodes( ElementNodes )
- n = GetElementNOFNodes()
- NodeIndexes => Element % NodeIndexes
- IF ( ANY( GM(GMPerm(NodeIndexes(:))) .GE. mskcrit ) ) CYCLE
- b = Basin(BasinPerm(e))
- IF (basinmax(b) < MAXVAL(distGL(distGLPerm(NodeIndexes(1:n))))) THEN
- basinmax(b)=MAXVAL(distGL(distGLPerm(NodeIndexes(1:n))))
- END IF
-
- IF (distmax < MAXVAL(distGL(distGLPerm(NodeIndexes(1:n))))) THEN
- distmax = MAXVAL(distGL(distGLPerm(NodeIndexes(1:n))))
- END IF
-END DO
-IF (Parallel) THEN
- CALL MPI_ALLREDUCE(distmax,max_Reduced,1,MPI_DOUBLE_PRECISION,MPI_MAX,ELMER_COMM_WORLD,ierr)
- distmax= max_Reduced
-END IF
-DO b=1,MaxBas
- IF (Parallel) THEN
- CALL MPI_ALLREDUCE(basinmax(b),basin_Reduced(b),1,MPI_DOUBLE_PRECISION,MPI_MAX,ELMER_COMM_WORLD,ierr)
- basinmax(b) = basin_Reduced(b)
- END IF
-END DO
-
-boxes = 1+NINT(SQRT(basinmax/distmax)*(boxmax-1))
-
-CALL INFO(TRIM(SolverName),'Boxes DONE', Level =5)
-
-!- Calculate total area of each box :
-! second loop on element
-DO e=1,Solver % NumberOfActiveElements
- Element => GetActiveElement(e)
- CALL GetElementNodes( ElementNodes )
- n = GetElementNOFNodes()
- NodeIndexes => Element % NodeIndexes
- VisitedNode(NodeIndexes(1:n))=VisitedNode(NodeIndexes(1:n))+1.0_dp
- ! leave the loop if grounded point in the element
- IF ( ANY( GM(GMPerm(NodeIndexes(:))) .GE. mskcrit ) ) CYCLE
- b= Basin(BasinPerm(e))
- nD=boxes(b)
- rr = SUM(distGL(distGLPerm(NodeIndexes(:))) &
- & / ( distGL(distGLPerm(NodeIndexes(:))) + distIF(distIFPerm(NodeIndexes(:))) )) &
- & / MAX(1,SIZE(distGL(distGLPerm(NodeIndexes(:)))))
- !localInteg = 0.0_dp
- localunity = 0.0_dp
- IntegStuff = GaussPoints( Element )
- DO t=1,IntegStuff % n
- U = IntegStuff % u(t)
- V = IntegStuff % v(t)
- W = IntegStuff % w(t)
- stat = ElementInfo(Element,ElementNodes,U,V,W,SqrtElementMetric, &
- Basis,dBasisdx )
- s = SqrtElementMetric * IntegStuff % s(t)
- localunity = localunity + s * SUM(Basis(1:n))
- END DO
- DO kk=1,nD
- IF ( rr .GE. 1.0-SQRT(1.0*(nD-kk+1)/nD) .AND. rr .LE. 1.0-SQRT(1.0*(nD-kk)/nD) ) THEN
- Abox(kk,b) = Abox(kk,b) + localunity
- ENDIF
- ENDDO
-END DO
-DO b=1,MaxBas
- nD=boxes(b)
- DO kk=1,nD
- IF (Parallel) THEN
- CALL MPI_ALLREDUCE(Abox(kk,b),Area_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
- Abox(kk,b) = Area_Reduced
- END IF
- ENDDO
-END DO
-
-! third loop on element
-Tbox(:,:)=0.d0 ; Sbox(:,:)=0.d0 ; qqq(:)=0.d0
-DO e=1,Solver % NumberOfActiveElements
- Element => GetActiveElement(e)
- CALL GetElementNodes( ElementNodes )
- n = GetElementNOFNodes()
- NodeIndexes => Element % NodeIndexes
- IF ( ANY( GM(GMPerm(NodeIndexes(:))) .GE. mskcrit ) ) CYCLE
- b= Basin(BasinPerm(e))
- nD=boxes(b)
- SELECT CASE(MeltVar % TYPE)
- CASE(Variable_on_elements)
- rr = SUM(distGL(distGLPerm(NodeIndexes(:))) &
- & / ( distGL(distGLPerm(NodeIndexes(:))) + distIF(distIFPerm(NodeIndexes(:))) )) &
- & / MAX(1,SIZE(distGL(distGLPerm(NodeIndexes(:)))))
- localunity = 0.0_dp
- IntegStuff = GaussPoints( Element )
- DO t=1,IntegStuff % n
- U = IntegStuff % u(t)
- V = IntegStuff % v(t)
- W = IntegStuff % w(t)
- stat = ElementInfo(Element,ElementNodes,U,V,W,SqrtElementMetric, &
- Basis,dBasisdx )
- s = SqrtElementMetric * IntegStuff % s(t)
- localunity = localunity + s * SUM(Basis(1:n))
- END DO
- IF ( rr .GE. 0.0 .AND. rr .LE. 1.0-SQRT(1.0*(nD-1)/nD) ) THEN
-
- zzz=SUM(Depth(DepthPerm(NodeIndexes(:))))/SIZE(NodeIndexes(:)) !mean depth of an element
- Tstar = lbd1*S0(b) + lbd2 + lbd3*zzz - T0(b) !NB: Tstar should be < 0
- g1 = gT * Abox(1,b)
- tmp1 = g1 / (CC*rhostar*(beta*S0(b)*meltfac-alpha))
- sn = (0.5*tmp1)**2 - tmp1*Tstar
- ! to avoid negative discriminent (no solution for x otherwise) :
- !if ( sn .lt. 0.d0 ) THEN
- ! xbox = 0.d0
- !else
- xbox = - 0.5*tmp1 + SQRT(sn) ! standard solution (Reese et al)
- !ENDif
- TT = T0(b) - xbox
- SS = S0(b) - xbox*S0(b)*meltfac
- Tbox(1,b) = Tbox(1,b) + TT * localunity
- Sbox(1,b) = Sbox(1,b) + SS * localunity
- qqq(b) = qqq(b) + CC*rhostar*(beta*(S0(b)-SS)-alpha*(T0(b)-TT)) * localunity
- Melt(MeltPerm(e)) = - gT * meltfac * ( lbd1*SS + lbd2 + lbd3*zzz - TT )
- Boxnumber(BPerm(e)) = 1
- print *,'z=',zzz,'Tstar=',Tstar,'g1=',g1,'tmp1=',tmp1,'sn=',sn,'xbox=',xbox,'Tbox=',Tbox(1,b),'Sbox=',Sbox(1,b),'qqq=',qqq(b),'Melt=', Melt(MeltPerm(e))
- CASE(Variable_on_nodes)
- DO nn=1,n
- zzz=Depth(DepthPerm(NodeIndexes(nn))) !depth of an node
- rr = distGL(distGLPerm((NodeIndexes(nn)))) / ( distGL(distGLPerm((NodeIndexes(nn)))) + distIF(distIFPerm((NodeIndexes(nn)))))
- IF ( rr .GE. 0.0 .AND. rr .LE. 1.0-SQRT(1.0*(nD-1)/nD) ) THEN
- Tstar = lbd1*S0(b) + lbd2 + lbd3*zzz - T0(b) !NB: Tstar should be < 0
- g1 = gT * Abox(1,b)
- tmp1 = g1 / (CC*rhostar*(beta*S0(b)*meltfac-alpha))
- sn = (0.5*tmp1)**2 - tmp1*Tstar
- ! to avoid negative discriminent (no solution for x otherwise) :
- !if ( sn .lt. 0.d0 ) THEN
- ! xbox = 0.d0
- !else
- xbox = - 0.5*tmp1 + SQRT(sn) ! standard solution (Reese et al)
- !ENDif
- TT = T0(b) - xbox
- SS = S0(b) - xbox*S0(b)*meltfac
- Tbox(1,b) = Tbox(1,b) + TT !* localunity
- Sbox(1,b) = Sbox(1,b) + SS !* localunity
- qqq(b) = qqq(b) + CC*rhostar*(beta*(S0(b)-SS)-alpha*(T0(b)-TT)) * localunity
-
- Melt(MeltPerm(NodeIndexes(nn))) = - gT * meltfac * ( lbd1*SS + lbd2 + lbd3*zzz - TT )
- Boxnumber(BPerm(NodeIndexes(nn))) = 1
- print *,'z=',zzz,'Tstar=',Tstar,'g1=',g1,'tmp1=',tmp1,'sn=',sn,'xbox=',xbox,'Tbox=',Tbox(1,b),'Sbox=',Sbox(1,b),'qqq=',qqq(b),'Melt=', Melt(MeltPerm(NodeIndexes(nn)))
- END IF
- END DO
- END SELECT
- END IF
-END DO
-DO b=1,MaxBas
- nD=boxes(b)
-
- IF (Parallel) THEN
- CALL MPI_ALLREDUCE(Tbox(1,b),Integ_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
- CALL MPI_ALLREDUCE(Sbox(1,b),Area_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
- CALL MPI_ALLREDUCE(qqq(b),qqq_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
- Tbox(1,b) = Integ_Reduced
- Sbox(1,b) = Area_Reduced
- qqq(b) = qqq_Reduced
- END IF
-ENDDO
-
-Tbox(1,1:MaxBas) = Tbox(1,1:MaxBas) / Abox(1,1:MaxBas)
-Sbox(1,1:MaxBas) = Sbox(1,1:MaxBas) / Abox(1,1:MaxBas)
-qqq(1:MaxBas) = qqq(1:MaxBas) / Abox(1,1:MaxBas)
-
-!- Temperature and salinity in possible other boxes :
-! 4 loops on the element
-DO kk=2,boxmax
-
- DO e=1,Solver % NumberOfActiveElements
- Element => GetActiveElement(e)
- CALL GetElementNodes( ElementNodes )
- n = GetElementNOFNodes()
- NodeIndexes => Element % NodeIndexes
- IF ( ANY( GM(GMPerm(NodeIndexes(:))) .GE. mskcrit ) ) CYCLE
- b= Basin(BasinPerm(e))
- nD=boxes(b)
- IF (kk .GT. nD) CYCLE
- SELECT CASE(MeltVar % TYPE)
- CASE(Variable_on_elements)
- rr = SUM(distGL(distGLPerm(NodeIndexes(:))) &
- & / ( distGL(distGLPerm(NodeIndexes(:))) + distIF(distIFPerm(NodeIndexes(:))))) &
- & / MAX(1,SIZE(distGL(distGLPerm(NodeIndexes(:)))))
- IF ( rr .GE. 1.0-SQRT(1.0*(nD-kk+1)/nD) .AND. rr .LE. 1.0-SQRT(1.0*(nD-kk)/nD) ) THEN
- localunity = 0.0_dp
- IntegStuff = GaussPoints( Element )
- DO t=1,IntegStuff % n
- U = IntegStuff % u(t)
- V = IntegStuff % v(t)
- W = IntegStuff % w(t)
- stat = ElementInfo(Element,ElementNodes,U,V,W,SqrtElementMetric, &
- Basis,dBasisdx )
- s = SqrtElementMetric * IntegStuff % s(t)
- localunity = localunity + s * SUM(Basis(1:n))
- END DO
- zzz = SUM(Depth(DepthPerm(NodeIndexes(:))))/SIZE(NodeIndexes(:)) !mean depth of an element
- Tstar = lbd1*Sbox(kk-1,b) + lbd2 + lbd3*zzz - Tbox(kk-1,b)
- g1 = gT * Abox(kk,b)
- g2 = g1 * meltfac
- xbox = - g1 * Tstar / ( qqq(b) + g1 - g2*lbd1*Sbox(kk-1,b) )
- TT = Tbox(kk-1,b) - xbox
- SS = Sbox(kk-1,b) - xbox*Sbox(kk-1,b)*meltfac
- Tbox(kk,b) = Tbox(kk,b) + TT * localunity
- Sbox(kk,b) = Sbox(kk,b) + SS * localunity
-
- Melt(MeltPerm(e)) = - gT * meltfac * ( lbd1*SS + lbd2 + lbd3*zzz - TT )
- Boxnumber(BPerm(e)) = kk
- print *,'kk=',kk,'z=',zzz,'Tstar=',Tstar,'g1=',g1,'g2=',g2,'xbox=',xbox,'Tbox=',Tbox(kk,b),'Sbox=',Sbox(kk,b),'Melt=', Melt(MeltPerm(e))
- CASE(Variable_on_nodes)
- DO nn=1,n
- zzz = Depth(DepthPerm(NodeIndexes(nn)))
- rr = distGL(distGLPerm((NodeIndexes(nn)))) / ( distGL(distGLPerm((NodeIndexes(nn)))) + distIF(distIFPerm((NodeIndexes(nn)))))
- IF ( rr .GE. 1.0-SQRT(1.0*(nD-kk+1)/nD) .AND. rr .LE. 1.0-SQRT(1.0*(nD-kk)/nD) ) THEN
- Tstar = lbd1*Sbox(kk-1,b) + lbd2 + lbd3*zzz - Tbox(kk-1,b)
- g1 = gT * Abox(kk,b)
- g2 = g1 * meltfac
- xbox = - g1 * Tstar / ( qqq(b) + g1 - g2*lbd1*Sbox(kk-1,b) )
- TT = Tbox(kk-1,b) - xbox
- SS = Sbox(kk-1,b) - xbox*Sbox(kk-1,b)*meltfac
- Tbox(kk,b) = Tbox(kk,b) + TT !* localunity
- Sbox(kk,b) = Sbox(kk,b) + SS !* localunity
- Melt(MeltPerm(NodeIndexes(nn))) = - gT * meltfac * ( lbd1*SS + lbd2 + lbd3*zzz - TT )
- Boxnumber(BPerm(NodeIndexes(nn))) = kk
- print *,'kk=',kk,'z=',zzz,'Tstar=',Tstar,'g1=',g1,'g2=',g2,'xbox=',xbox,'Tbox=',Tbox(kk,b),'Sbox=',Sbox(kk,b),'Melt=', Melt(MeltPerm(NodeIndexes(nn)))
- END IF
- END DO
- END SELECT
- ENDif
- END DO
- DO b=1,MaxBas
- nD=boxes(b)
- IF (Parallel) THEN
- CALL MPI_ALLREDUCE(Tbox(kk,b),Integ_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
- CALL MPI_ALLREDUCE(Sbox(kk,b),Area_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
- Tbox(kk,b) = Integ_Reduced
- Sbox(kk,b) = Area_Reduced
- END IF
- ENDDO
- Tbox(kk,1:MaxBas) = Tbox(kk,1:MaxBas) / Abox(kk,1:MaxBas)
- Sbox(kk,1:MaxBas) = Sbox(kk,1:MaxBas) / Abox(kk,1:MaxBas)
-END DO
+ ! first loop on element to determine the number of boxes per basins
+ DO e=1,Solver % NumberOfActiveElements
+ Element => GetActiveElement(e)
+ CALL GetElementNodes( ElementNodes )
+ n = GetElementNOFNodes()
+ NodeIndexes => Element % NodeIndexes
+ IF ( ANY( GM(GMPerm(NodeIndexes(:))) .GE. mskcrit ) ) CYCLE
+ b = Basin(BasinPerm(e))
+ IF (basinmax(b) < MAXVAL(distGL(distGLPerm(NodeIndexes(1:n))))) THEN
+ basinmax(b)=MAXVAL(distGL(distGLPerm(NodeIndexes(1:n))))
+ END IF
+
+ IF (distmax < MAXVAL(distGL(distGLPerm(NodeIndexes(1:n))))) THEN
+ distmax = MAXVAL(distGL(distGLPerm(NodeIndexes(1:n))))
+ END IF
+ END DO
+ IF (Parallel) THEN
+ CALL MPI_ALLREDUCE(distmax,max_Reduced,1,MPI_DOUBLE_PRECISION,MPI_MAX,ELMER_COMM_WORLD,ierr)
+ distmax= max_Reduced
+ END IF
+ DO b=1,MaxBas
+ IF (Parallel) THEN
+ CALL MPI_ALLREDUCE(basinmax(b),basin_Reduced(b),1,MPI_DOUBLE_PRECISION,MPI_MAX,ELMER_COMM_WORLD,ierr)
+ basinmax(b) = basin_Reduced(b)
+ END IF
+ END DO
+
+ boxes = 1+NINT(SQRT(basinmax/distmax)*(boxmax-1))
+
+ CALL INFO(TRIM(SolverName),'Boxes DONE', Level =5)
+
+ !- Calculate total area of each box :
+ ! second loop on element
+ DO e=1,Solver % NumberOfActiveElements
+ Element => GetActiveElement(e)
+ CALL GetElementNodes( ElementNodes )
+ n = GetElementNOFNodes()
+ NodeIndexes => Element % NodeIndexes
+ VisitedNode(NodeIndexes(1:n))=VisitedNode(NodeIndexes(1:n))+1.0_dp
+ ! leave the loop if grounded point in the element
+ IF ( ANY( GM(GMPerm(NodeIndexes(:))) .GE. mskcrit ) ) CYCLE
+ b= Basin(BasinPerm(e))
+ nD=boxes(b)
+ rr = SUM(distGL(distGLPerm(NodeIndexes(:))) &
+ & / ( distGL(distGLPerm(NodeIndexes(:))) + distIF(distIFPerm(NodeIndexes(:))) )) &
+ & / MAX(1,SIZE(distGL(distGLPerm(NodeIndexes(:)))))
+ !localInteg = 0.0_dp
+ localunity = 0.0_dp
+ IntegStuff = GaussPoints( Element )
+ DO t=1,IntegStuff % n
+ U = IntegStuff % u(t)
+ V = IntegStuff % v(t)
+ W = IntegStuff % w(t)
+ stat = ElementInfo(Element,ElementNodes,U,V,W,SqrtElementMetric, &
+ Basis,dBasisdx )
+ s = SqrtElementMetric * IntegStuff % s(t)
+ localunity = localunity + s * SUM(Basis(1:n))
+ END DO
+ DO kk=1,nD
+ IF ( rr .GE. 1.0-SQRT(1.0*(nD-kk+1)/nD) .AND. rr .LE. 1.0-SQRT(1.0*(nD-kk)/nD) ) THEN
+ Abox(kk,b) = Abox(kk,b) + localunity
+ ENDIF
+ ENDDO
+ END DO
+ DO b=1,MaxBas
+ nD=boxes(b)
+ DO kk=1,nD
+ IF (Parallel) THEN
+ CALL MPI_ALLREDUCE(Abox(kk,b),Area_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
+ Abox(kk,b) = Area_Reduced
+ END IF
+ ENDDO
+ END DO
+
+ ! third loop on element
+ Tbox(:,:)=0.d0 ; Sbox(:,:)=0.d0 ; qqq(:)=0.d0
+ DO e=1,Solver % NumberOfActiveElements
+ Element => GetActiveElement(e)
+ CALL GetElementNodes( ElementNodes )
+ n = GetElementNOFNodes()
+ NodeIndexes => Element % NodeIndexes
+ IF ( ANY( GM(GMPerm(NodeIndexes(:))) .GE. mskcrit ) ) CYCLE
+ b= Basin(BasinPerm(e))
+ nD=boxes(b)
+ SELECT CASE(MeltVar % TYPE)
+ CASE(Variable_on_elements)
+ rr = SUM(distGL(distGLPerm(NodeIndexes(:))) &
+ & / ( distGL(distGLPerm(NodeIndexes(:))) + distIF(distIFPerm(NodeIndexes(:))) )) &
+ & / MAX(1,SIZE(distGL(distGLPerm(NodeIndexes(:)))))
+ localunity = 0.0_dp
+ IntegStuff = GaussPoints( Element )
+ DO t=1,IntegStuff % n
+ U = IntegStuff % u(t)
+ V = IntegStuff % v(t)
+ W = IntegStuff % w(t)
+ stat = ElementInfo(Element,ElementNodes,U,V,W,SqrtElementMetric, &
+ Basis,dBasisdx )
+ s = SqrtElementMetric * IntegStuff % s(t)
+ localunity = localunity + s * SUM(Basis(1:n))
+ END DO
+ IF ( rr .GE. 0.0 .AND. rr .LE. 1.0-SQRT(1.0*(nD-1)/nD) ) THEN
+
+ zzz=SUM(Depth(DepthPerm(NodeIndexes(:))))/SIZE(NodeIndexes(:)) !mean depth of an element
+ Tstar = lbd1*S0(b) + lbd2 + lbd3*zzz - T0(b) !NB: Tstar should be < 0
+ g1 = gT * Abox(1,b)
+ tmp1 = g1 / (CC*rhostar*(beta*S0(b)*meltfac-alpha))
+ sn = (0.5*tmp1)**2 - tmp1*Tstar
+ ! to avoid negative discriminent (no solution for x otherwise) :
+ !if ( sn .lt. 0.d0 ) THEN
+ ! xbox = 0.d0
+ !else
+ xbox = - 0.5*tmp1 + SQRT(sn) ! standard solution (Reese et al)
+ !ENDif
+ TT = T0(b) - xbox
+ SS = S0(b) - xbox*S0(b)*meltfac
+ Tbox(1,b) = Tbox(1,b) + TT * localunity
+ Sbox(1,b) = Sbox(1,b) + SS * localunity
+ qqq(b) = qqq(b) + CC*rhostar*(beta*(S0(b)-SS)-alpha*(T0(b)-TT)) * localunity
+ Melt(MeltPerm(e)) = - gT * meltfac * ( lbd1*SS + lbd2 + lbd3*zzz - TT )
+ Boxnumber(BPerm(e)) = 1
+ print *,'z=',zzz,'Tstar=',Tstar,'g1=',g1,'tmp1=',tmp1,'sn=',sn,'xbox=',xbox,'Tbox=',Tbox(1,b),'Sbox=',Sbox(1,b),'qqq=',qqq(b),'Melt=', Melt(MeltPerm(e))
+ END IF
+ CASE(Variable_on_nodes)
+ DO nn=1,n
+ zzz=Depth(DepthPerm(NodeIndexes(nn))) !depth of an node
+ rr = distGL(distGLPerm((NodeIndexes(nn)))) / ( distGL(distGLPerm((NodeIndexes(nn)))) + distIF(distIFPerm((NodeIndexes(nn)))))
+ IF ( rr .GE. 0.0 .AND. rr .LE. 1.0-SQRT(1.0*(nD-1)/nD) ) THEN
+ Tstar = lbd1*S0(b) + lbd2 + lbd3*zzz - T0(b) !NB: Tstar should be < 0
+ g1 = gT * Abox(1,b)
+ tmp1 = g1 / (CC*rhostar*(beta*S0(b)*meltfac-alpha))
+ sn = (0.5*tmp1)**2 - tmp1*Tstar
+ ! to avoid negative discriminent (no solution for x otherwise) :
+ !if ( sn .lt. 0.d0 ) THEN
+ ! xbox = 0.d0
+ !else
+ xbox = - 0.5*tmp1 + SQRT(sn) ! standard solution (Reese et al)
+ !ENDif
+ TT = T0(b) - xbox
+ SS = S0(b) - xbox*S0(b)*meltfac
+ Tbox(1,b) = Tbox(1,b) + TT !* localunity
+ Sbox(1,b) = Sbox(1,b) + SS !* localunity
+ qqq(b) = qqq(b) + CC*rhostar*(beta*(S0(b)-SS)-alpha*(T0(b)-TT))! * localunity
+
+ Melt(MeltPerm(NodeIndexes(nn))) = - gT * meltfac * ( lbd1*SS + lbd2 + lbd3*zzz - TT )
+ Boxnumber(BPerm(NodeIndexes(nn))) = 1
+ print *,'z=',zzz,'Tstar=',Tstar,'g1=',g1,'tmp1=',tmp1,'sn=',sn,'xbox=',xbox,'Tbox=',Tbox(1,b),'Sbox=',Sbox(1,b),'qqq=',qqq(b),'Melt=', Melt(MeltPerm(NodeIndexes(nn)))
+ END IF
+ END DO
+ END SELECT
+ END DO
+ DO b=1,MaxBas
+ nD=boxes(b)
+
+ IF (Parallel) THEN
+ CALL MPI_ALLREDUCE(Tbox(1,b),Integ_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
+ CALL MPI_ALLREDUCE(Sbox(1,b),Area_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
+ CALL MPI_ALLREDUCE(qqq(b),qqq_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
+ Tbox(1,b) = Integ_Reduced
+ Sbox(1,b) = Area_Reduced
+ qqq(b) = qqq_Reduced
+ END IF
+ ENDDO
+
+ Tbox(1,1:MaxBas) = Tbox(1,1:MaxBas) / Abox(1,1:MaxBas)
+ Sbox(1,1:MaxBas) = Sbox(1,1:MaxBas) / Abox(1,1:MaxBas)
+ qqq(1:MaxBas) = qqq(1:MaxBas) / Abox(1,1:MaxBas)
+
+ !- Temperature and salinity in possible other boxes :
+ ! 4 loops on the element
+ DO kk=2,boxmax
+
+ DO e=1,Solver % NumberOfActiveElements
+ Element => GetActiveElement(e)
+ CALL GetElementNodes( ElementNodes )
+ n = GetElementNOFNodes()
+ NodeIndexes => Element % NodeIndexes
+ IF ( ANY( GM(GMPerm(NodeIndexes(:))) .GE. mskcrit ) ) CYCLE
+ b= Basin(BasinPerm(e))
+ nD=boxes(b)
+ IF (kk .GT. nD) CYCLE
+ SELECT CASE(MeltVar % TYPE)
+ CASE(Variable_on_elements)
+ rr = SUM(distGL(distGLPerm(NodeIndexes(:))) &
+ & / ( distGL(distGLPerm(NodeIndexes(:))) + distIF(distIFPerm(NodeIndexes(:))))) &
+ & / MAX(1,SIZE(distGL(distGLPerm(NodeIndexes(:)))))
+ IF ( rr .GE. 1.0-SQRT(1.0*(nD-kk+1)/nD) .AND. rr .LE. 1.0-SQRT(1.0*(nD-kk)/nD) ) THEN
+ localunity = 0.0_dp
+ IntegStuff = GaussPoints( Element )
+ DO t=1,IntegStuff % n
+ U = IntegStuff % u(t)
+ V = IntegStuff % v(t)
+ W = IntegStuff % w(t)
+ stat = ElementInfo(Element,ElementNodes,U,V,W,SqrtElementMetric, &
+ Basis,dBasisdx )
+ s = SqrtElementMetric * IntegStuff % s(t)
+ localunity = localunity + s * SUM(Basis(1:n))
+ END DO
+ zzz = SUM(Depth(DepthPerm(NodeIndexes(:))))/SIZE(NodeIndexes(:)) !mean depth of an element
+ Tstar = lbd1*Sbox(kk-1,b) + lbd2 + lbd3*zzz - Tbox(kk-1,b)
+ g1 = gT * Abox(kk,b)
+ g2 = g1 * meltfac
+ xbox = - g1 * Tstar / ( qqq(b) + g1 - g2*lbd1*Sbox(kk-1,b) )
+ TT = Tbox(kk-1,b) - xbox
+ SS = Sbox(kk-1,b) - xbox*Sbox(kk-1,b)*meltfac
+ Tbox(kk,b) = Tbox(kk,b) + TT * localunity
+ Sbox(kk,b) = Sbox(kk,b) + SS * localunity
+
+ Melt(MeltPerm(e)) = - gT * meltfac * ( lbd1*SS + lbd2 + lbd3*zzz - TT )
+ Boxnumber(BPerm(e)) = kk
+ print *,'kk=',kk,'z=',zzz,'Tstar=',Tstar,'g1=',g1,'g2=',g2,'xbox=',xbox,'Tbox=',Tbox(kk,b),'Sbox=',Sbox(kk,b),'Melt=', Melt(MeltPerm(e))
+ END IF
+
+ CASE(Variable_on_nodes)
+ DO nn=1,n
+ zzz = Depth(DepthPerm(NodeIndexes(nn)))
+ rr = distGL(distGLPerm((NodeIndexes(nn)))) / ( distGL(distGLPerm((NodeIndexes(nn)))) + distIF(distIFPerm((NodeIndexes(nn)))))
+ IF ( rr .GE. 1.0-SQRT(1.0*(nD-kk+1)/nD) .AND. rr .LE. 1.0-SQRT(1.0*(nD-kk)/nD) ) THEN
+ Tstar = lbd1*Sbox(kk-1,b) + lbd2 + lbd3*zzz - Tbox(kk-1,b)
+ g1 = gT * Abox(kk,b)
+ g2 = g1 * meltfac
+ xbox = - g1 * Tstar / ( qqq(b) + g1 - g2*lbd1*Sbox(kk-1,b) )
+ TT = Tbox(kk-1,b) - xbox
+ SS = Sbox(kk-1,b) - xbox*Sbox(kk-1,b)*meltfac
+ Tbox(kk,b) = Tbox(kk,b) + TT !* localunity
+ Sbox(kk,b) = Sbox(kk,b) + SS !* localunity
+ Melt(MeltPerm(NodeIndexes(nn))) = - gT * meltfac * ( lbd1*SS + lbd2 + lbd3*zzz - TT )
+ Boxnumber(BPerm(NodeIndexes(nn))) = kk
+ print *,'kk=',kk,'z=',zzz,'Tstar=',Tstar,'g1=',g1,'g2=',g2,'xbox=',xbox,'Tbox=',Tbox(kk,b),'Sbox=',Sbox(kk,b),'Melt=', Melt(MeltPerm(NodeIndexes(nn)))
+ END IF
+ END DO
+ END SELECT
+ END DO
+ DO b=1,MaxBas
+ nD=boxes(b)
+ IF (Parallel) THEN
+ CALL MPI_ALLREDUCE(Tbox(kk,b),Integ_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
+ CALL MPI_ALLREDUCE(Sbox(kk,b),Area_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
+ Tbox(kk,b) = Integ_Reduced
+ Sbox(kk,b) = Area_Reduced
+ END IF
+ ENDDO
+ Tbox(kk,1:MaxBas) = Tbox(kk,1:MaxBas) / Abox(kk,1:MaxBas)
+ Sbox(kk,1:MaxBas) = Sbox(kk,1:MaxBas) / Abox(kk,1:MaxBas)
+ END DO
!!! last loop on the element
!!! USEFULL or USELESS to compute Integrale of the melt?????
-Integ = 0.0_dp
-DO e=1,Solver % NumberOfActiveElements
- Element => GetActiveElement(e)
- CALL GetElementNodes( ElementNodes )
- n = GetElementNOFNodes(Element)
- NodeIndexes => Element % NodeIndexes
- VisitedNode(NodeIndexes(1:n))=VisitedNode(NodeIndexes(1:n))+1.0_dp
- localInteg = 0.0_dp
- IntegStuff = GaussPoints( Element )
- DO t=1,IntegStuff % n
- U = IntegStuff % u(t)
- V = IntegStuff % v(t)
- W = IntegStuff % w(t)
- stat = ElementInfo(Element,ElementNodes,U,V,W,SqrtElementMetric, &
- Basis,dBasisdx )
- s = SqrtElementMetric * IntegStuff % s(t)
- SELECT CASE(MeltVar % TYPE)
- CASE(Variable_on_elements)
- localInteg = localInteg + s * SUM(Basis(:) * Melt(MeltPerm(e)))
- CASE(Variable_on_nodes)
- localInteg = localInteg + s * SUM(Basis(:) * Melt(MeltPerm(NodeIndexes(:))))
- END SELECT
-
-
- END DO
- Integ = Integ + localInteg
-ENDDO
-
-IF (Parallel) THEN
- CALL MPI_ALLREDUCE(Integ,Integ_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
- IF (Solver % Matrix % ParMatrix % ParEnv % MyPE == 0) THEN
- WRITE(meltValue,'(F20.2)') Integ_Reduced
- Message='TOTAL_MELT_RATE: '//meltValue
- CALL INFO(SolverName,Message,Level=1)
- END IF
-ELSE
- Integ_Reduced = Integ
-ENDIF
-
-DEALLOCATE( Zbox, Abox, Tbox, Sbox, Mbox, T0, S0)
-
-DEALLOCATE(basin_Reduced,basinmax,boxes)
-DEALLOCATE(VisitedNode, Basis, dBasisdx)
+ Integ = 0.0_dp
+ DO e=1,Solver % NumberOfActiveElements
+ Element => GetActiveElement(e)
+ CALL GetElementNodes( ElementNodes )
+ n = GetElementNOFNodes(Element)
+ NodeIndexes => Element % NodeIndexes
+ VisitedNode(NodeIndexes(1:n))=VisitedNode(NodeIndexes(1:n))+1.0_dp
+ localInteg = 0.0_dp
+ IntegStuff = GaussPoints( Element )
+ DO t=1,IntegStuff % n
+ U = IntegStuff % u(t)
+ V = IntegStuff % v(t)
+ W = IntegStuff % w(t)
+ stat = ElementInfo(Element,ElementNodes,U,V,W,SqrtElementMetric, &
+ Basis,dBasisdx )
+ s = SqrtElementMetric * IntegStuff % s(t)
+ SELECT CASE(MeltVar % TYPE)
+ CASE(Variable_on_elements)
+ localInteg = localInteg + s * SUM(Basis(:) * Melt(MeltPerm(e)))
+ CASE(Variable_on_nodes)
+ localInteg = localInteg + s * SUM(Basis(:) * Melt(MeltPerm(NodeIndexes(:))))
+ END SELECT
+
+
+ END DO
+ Integ = Integ + localInteg
+ ENDDO
+
+ IF (Parallel) THEN
+ CALL MPI_ALLREDUCE(Integ,Integ_Reduced,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
+ IF (Solver % Matrix % ParMatrix % ParEnv % MyPE == 0) THEN
+ WRITE(meltValue,'(F20.2)') Integ_Reduced
+ Message='TOTAL_MELT_RATE: '//meltValue
+ CALL INFO(SolverName,Message,Level=1)
+ END IF
+ ELSE
+ Integ_Reduced = Integ
+ ENDIF
+
+ DEALLOCATE( Zbox, Abox, Tbox, Sbox, Mbox, T0, S0)
+
+ DEALLOCATE(basin_Reduced,basinmax,boxes)
+ DEALLOCATE(VisitedNode, Basis, dBasisdx)
END SUBROUTINE PICO
diff --git a/src/boxmodel_solver.md b/src/boxmodel_solver.md
index 22265be..08be4fb 100644
--- a/src/boxmodel_solver.md
+++ b/src/boxmodel_solver.md
@@ -1,5 +1,7 @@
# Boxes model Solver
+
+
## General information
- **Solver Fortran File:** BoxModel.F90
@@ -12,7 +14,7 @@
## General Description
-This solver is devellop to compute the basal melt by using the Box Model (Reese, 2018). It will be used with the nearestpoint solver to have the basin variable. To compute the distance variable, the disance solver and FrontThickness_mask solver are used.
+This solver is develop to compute the basal melt by using the Box Model (Reese, 2018). It will be used with the nearestpoint solver to have the basin variable. To compute the distance variable, the distance solver and FrontThickness_mask solver are used.
In the sif, many constant have to be define. This constant aren’t specific to the box model
@@ -65,7 +67,7 @@ Box boundary :
-- Then the solver compute the **total area of each boxe**
+- Then the solver compute the **total area of each box**
- The goal is to solve the melt equation :
@@ -186,3 +188,293 @@ End
```
+## Algorithme
+
+
+
+### Définition du nombre de boite
+
+basinmax(:) = 0 # vecteur de longueur du nombre de basin
+
+distmax = 0
+
+Abox(:,:)=0 #matrice nombre de bassin x nombre de boite
+
+melt(:)=0.0dp
+
+Boucle sur les élément :
+
+ Si un node de l’élément GM >= mskcrit alors passer à l’itération suivante
+
+ b <- valeur bassin de l’élément
+
+ si basinmax(b) < max(distGL) alors
+
+ basinmax(b) = max(distGL)
+
+ FIN si
+
+ si distmax < max(distGL)
+
+ distmax = max(distGL)
+
+ FIN si
+
+FIN boucle
+
+$$boxes = 1 + rd(\sqrt{basinlax/distmax}*(n_{max}-1))$$ # nombre de boite par bassin / vecteur de longueur du nombre de basin
+
+### Calcul de l’aire des boites
+
+Abox (:,:) =0.0
+
+Boucle sur les éléments e= 1 , élément :
+
+ node = noeuds de e
+
+ Si un node de l’élément GM(node) >= mskcrit alors passer à l’itération suivante
+
+ b <- valeur bassin(e) de l’élément
+
+ nD = boxes (b) # définition du nombre de boite pour le bassin b
+
+ $$rr = \sum\frac{distGL(node)}{distGL(node) + distIF(node)} / len(distGL(node))$$# calcul de la distance relative de l’élément au front
+
+ localunity = 0.0
+
+ Boucle sur les points d’intégration de l’élément :
+
+ définition du poids d’intégration **s**
+
+ localunity = localunity + s*SUM(BAsis(1:n)) #calcul de l’aire de l’élément
+
+ FIN boucle
+
+ Boucle sur le nombre de boite du bassin kk = 1 à nD :
+
+ Si $$1.0-\sqrt{1.0*(n_{max}-kk+1)/n_{max}} < rr < 1.0-\sqrt{1.0*(n_{max}-kk)/n_{max}} $$ alors
+
+ Abox(kk,b) = Abox(kk,b) +localunity #aire des boites par bassins
+
+ FIN si
+
+ FIN Boucle
+
+FIN Boucle
+
+### Calcul fonte dans la première boite
+
+Tbox(:,:) =0.0
+
+Sbox(:,:)=0.0
+
+qqq(:)=0.0
+
+Boucle sur les éléments e= 1 , element:
+
+ node = noeuds de e
+
+ Si un node de l’élément GM(node) >= mskcrit alors passer à l’itération suivante
+
+ b <- valeur bassin(e) de l’élément
+
+ nD = boxes (b) # définition du nombre de boite pour le bassin b
+
+ **CAS 1** sur les éléments :
+
+ $$rr = \sum\frac{distGL(node)}{distGL(node) + distIF(node)} / len(distGL(node))$$ # calcul de la distance relative de l’élément au front
+
+ localunity = 0.0
+
+ Boucle sur les points d’intégration de l’élément n:
+
+ définition du poids d’intégration s
+
+ localunity = localunity + s*SUM(BAsis(1:n)) #calcul de l’aire de l’élément
+
+ FIN boucle
+
+ Si $$0 < rr < 1.0-\sqrt{1.0*(n_{max}-1)/n_{max}} $$ alors
+
+ $$z_b$$ profondeur de l’élément
+
+ $$T_*=aS_0(bb)+b+c*z_b-T_0(bb)$$ # calcul T* difference entre la température de fonte et la temperature de la cavité
+
+ $$g1 = \gamma^*_T * A_{box}(1,bb)$$ # calcul du flux d’échange de la boite 1
+
+ $$tmp1 =g1/(C*\rho_**(\beta*S_0*meltfac-\alpha)) $$
+
+ $$xbox=0.5*tmp1+ \sqrt{(0.5*tmp1)^2-tmp1*T_*}$$
+
+ $$TT=T_0(bb)-xbox$$ # Calcul de l’évolution de la temperature entre la boite 0 et 1 pour un élément
+
+ $$SS=S_0(bb)-xbox*S_0(bb)*meltfac $$ # Calcul de l’évolution de la salinité entre la boite 0 et 1 pour un élément
+
+ $$T(1,bb)=T(1,bb) + TT * localunity$$ # évolution de la température pondéré par l’aire de l’éléments
+
+ $$S(1,bb)=S(1,bb) + SS * localunity$$ # évolution de la salinité pondéré par l’aire de l’éléemnts
+
+ $$qqq(bb) = qqq(bb) + C \rho_*(\beta(S_{0}(bb)-SS)-\alpha(T_{0}(bb)-TT)) * localunity$$ #overturning flux moyen de boite 1
+
+ $$Melt(e) = -\gamma^*_T*meltfac* (aSS+b+c*z_b-TT)$$
+
+ FIN Si
+
+ **CAS 2** sur les noeuds :
+
+ Boucle sur node = nn :
+
+ $z_b$ profondeur du noeud
+
+ calcul de $$rr = \frac{distGL(nn)}{distGL(nn) + distIF(nn)} $$
+
+ Si $$0 < rr < 1.0-\sqrt{1.0*(n_{max}-1)/n_{max}} $$ alors
+
+ $$T_*=aS_0+b+c*z_b-T_0$$ # calcul T*
+
+ $$g1 = \gamma^*_T * A_{box}(1,bb)$$
+
+ $$tmp1 =g1/(C*\rho_**(\beta*S_0(bb)*meltfac-\alpha)) $$
+
+ $$xbox=0.5*tmp1+ \sqrt{(0.5*tmp1)^2-tmp1*T_*}$$
+
+ $$TT=T_0(bb)-xbox$$
+
+ $$SS=S_0(bb)-xbox*S_0(bb)*meltfac$$
+
+ $$T(1,bb)=T(1,bb)+TT$$
+
+ $$S(1,bb)=S(1,bb)+SS$$
+
+ $$qqq(bb) = qqq(bb) + C \rho_*(\beta(S_{0}(bb)-SS)-\alpha(T_{0}(bb)-TT)$$
+
+ $$Melt(nn) = -\gamma^*_T*meltfac* (aSS+b+c*z_b-TT)$$
+
+
+
+ FIN Si
+
+ FIN Boucle
+
+Fin Boucle
+
+Boucle sur les basins bb
+
+ SI parrallel = True Alors
+
+ Faire mpi_somme de Tbox(1,bb)
+
+ Faire mpi_somme de Abox(1,bb)
+
+ Faire somme de qqq(bb)
+
+ FIN si
+
+FIN boucle
+
+### Calcul de la fonte pour les boite 2 à $$n_{max}$$
+
+Boucle sur les boites kk = 2 à $$n_{max}$$
+
+ Boucle sur e= 1 , element:
+
+ node = noeud de e
+
+ Si un node de l’élément GM(node) >= mskcrit alors passer à l’itération suivante
+
+ bb <- valeur bassin(e) de l’élément
+
+ nD = boxes (b) # définition du nombre de boite pour le bassin b
+
+ Si kk > nD passer à l’itération suivante
+
+ **CAS 1** sur les éléments :
+
+ nD = boxes (b) # définition du nombre de boite pour le bassin b
+
+ $$rr = \sum\frac{distGL(node)}{distGL(node) + distIF(node)} / len(distGL(node))$$ # calcul de la distance relative de l’élément au front
+
+ localunity = 0.0
+
+ Boucle sur les points d’intégration de l’élément :
+
+ définition du poids d’intégration s
+
+ localunity = localunity + s*SUM(BAsis(1:n)) #calcul de l’aire de l’élément
+
+ FIN boucle
+
+ Si $$1.0-\sqrt{1.0*(n_{max}-kk+1)/n_{max}} < rr < 1.0-\sqrt{1.0*(n_{max}-kk)/n_{max}} $$ alors
+
+ $$z_b$$ profondeur de l’élément
+
+ $$T_*=aS(kk-1,bb)+b+c*z_b-T_{k-1}$$
+
+ $$g1 = gT * A_{box}(kk,bb)$$
+
+ $$g2 = g1 * meltfac$$
+
+ $$xbox=- g1 * T_* / ( q + g1 - g2*a*S_{k-1} )$$
+
+ $$TT = T(kk-1,bb)-xbox$$
+
+ $$SS = S(kk-1,bb)-xbox*S(kk-1,bb)*meltfac$$
+
+ $$T(kk,bb)=T(kk,bb)+TT * localunity$$
+
+ $$S(kk,bb)=S(kk,bb)+SS*localunity$$
+
+ $$Melt(e) = \gamma^*_T*meltfac* (aS_k+b-c*z_b-T_k)$$
+
+ FIN Si
+
+ **CAS 2** sur les nœuds :
+
+ Boucle sur les nœud de l’élément :
+
+ $z_b$ profondeur du nœud
+
+ calcul de $$rr = \frac{distGL(node)}{distGL(node) + distIF(node)} $$
+
+ Si $$ 1.0-\sqrt{1.0*(n_{max}-kk+1)/n_{max}}< rr < 1.0-\sqrt{1.0*(n_{max}-kk)/n_{max}} $$ alors
+
+ $$T_*=aS_{k-1}+b+c*z_b-T_{k-1}$$
+
+ $$g1 = gT * A_{box}(kk,bb)$$
+
+ $$g2 = g1 * meltfac$$
+
+ $$xbox=- g1 * T_* / ( q + g1 - g2*a*S_{k-1} )$$
+
+ $$TT = T(kk-1,bb)-xbox$$
+
+ $$SS = S(kk-1,bb)-xbox*S(kk-1,bb)*meltfac$$
+
+ $$T(kk,bb)=T(kk,bb)+TT$$
+
+ $$S(kk,bb)=S(kk,bb)+SS$$
+
+ $$Melt(node) = \gamma^*_T*meltfac* (aSS+b-c*z_b-TT)$$
+
+ FIN Si
+
+ FIN Boucle
+
+ Fin Boucle
+
+ Boucle sur les basins bb
+
+ SI parrallel = True Alors
+
+ Faire mpi_somme de Tbox(kk,bb)
+
+ Faire mpi_somme de Abox(kk,bb)
+
+ FIN si
+
+ FIN boucle
+
+FIN Boucle
+
+
+
--
GitLab