Skip to content

Commit

Permalink
Enable pickup for pH3D for use with calcite saturation calculation
Browse files Browse the repository at this point in the history
  • Loading branch information
seamanticscience committed Apr 5, 2023
1 parent b9d4fbd commit a7cfa71
Show file tree
Hide file tree
Showing 6 changed files with 582 additions and 93 deletions.
3 changes: 2 additions & 1 deletion pkg/dic/DIC_VARS.h
Original file line number Diff line number Diff line change
Expand Up @@ -68,8 +68,9 @@ C silicaDeep :: 3D-field of silicate concentration for pH and
C carbonate calculations. Read in from file (mol/m3).
C omegaC :: Local saturation state with respect to calcite
COMMON /DIC_CALCITE_SAT_FIELDS/
& silicaDeep, omegaC
& silicaDeep, pH3D, omegaC
_RL silicaDeep(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL pH3D (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL omegaC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
#endif

Expand Down
47 changes: 25 additions & 22 deletions pkg/dic/calcite_saturation.F
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ SUBROUTINE CALCITE_SATURATION( PTR_DIC, PTR_ALK, PTR_PO4,
locSalt(i,j) = salt (i,j,k,bi,bj)
carbonate(i,j) = 0. _d 0
pCO2local(i,j) = 0. _d 0
pHlocal(i,j) = 0. _d 0
pHlocal(i,j) = pH3D(i,j,k,bi,bj)
ENDDO
ENDDO
#ifdef CARBONCHEM_SOLVESAPHE
Expand Down Expand Up @@ -94,7 +94,6 @@ SUBROUTINE CALCITE_SATURATION( PTR_DIC, PTR_ALK, PTR_PO4,
debugPrt = debugMode
DO j=jmin,jmax
DO i=imin,imax

IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
#ifdef CARBONCHEM_SOLVESAPHE
C calcium concentration already calculated in S/R DIC_COEFFS_SURF
Expand All @@ -113,15 +112,16 @@ SUBROUTINE CALCITE_SATURATION( PTR_DIC, PTR_ALK, PTR_PO4,
#ifdef CARBONCHEM_SOLVESAPHE
IF ( selectPHsolver.GT.0 ) THEN
C Use Munhoven (2013) Solvesaphe routine to calculate pH and pCO2
#ifdef ALLOW_DEBUG
IF (debugPrt) CALL DEBUG_CALL('AHINI_FOR_AT',myThid)
#endif
C call AHINI_FOR_AT to get better initial guess of pH
CALL AHINI_FOR_AT(alklocal*permil,
I diclocal*permil,
I bt(i,j,bi,bj),
O pHlocal(i,j),
I i,j,k,bi,bj,myIter,myThid )
C#ifdef ALLOW_DEBUG
C IF (debugPrt) CALL DEBUG_CALL('AHINI_FOR_AT',myThid)
C#endif
C This is moved to dic_ini_forcing
CC call AHINI_FOR_AT to get better initial guess of pH
C CALL AHINI_FOR_AT(alklocal*permil,
C I diclocal*permil,
C I bt(i,j,bi,bj),
C O pHlocal(i,j),
C I i,j,k,bi,bj,myIter,myThid )

#ifdef ALLOW_DEBUG
IF (debugPrt) CALL DEBUG_CALL('CALC_PCO2_SOLVESAPHE',myThid)
Expand All @@ -141,13 +141,13 @@ SUBROUTINE CALCITE_SATURATION( PTR_DIC, PTR_ALK, PTR_PO4,
#ifdef ALLOW_DEBUG
IF (debugPrt) CALL DEBUG_CALL('CALC_PCO2_APPROX',myThid)
#endif
pHlocal(i,j) = 7.9 _d 0

C Since we start cold with an arbitrary pH, we must iterate pH solver
C to ensure accurate estimate of co3 at depth
C Because we may call this for deep ocean infrequently can afford to make
C several iterations

C This is moved to dic_ini_forcing, so pH is no longer a cold start
CC pHlocal(i,j) = 7.9 _d 0
CC Since we start cold with an arbitrary pH, we must iterate pH solver
CC to ensure accurate estimate of co3 at depth
CC Because we may call this for deep ocean infrequently can afford to make
CC several iterations
C But we could still iterate, especially if calcOmegaCalciteFreq is more than a few timesteps
DO CO3iter = 1, nIterCO3
CALL CALC_PCO2_APPROX(
I locTemp(i,j), locSalt(i,j),
Expand All @@ -158,9 +158,9 @@ SUBROUTINE CALCITE_SATURATION( PTR_DIC, PTR_ALK, PTR_PO4,
I aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
I aksi(i,j,bi,bj),akf(i,j,bi,bj),
I ak0(i,j,bi,bj), fugf(i,j,bi,bj), ff(i,j,bi,bj),
I bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
I bt(i,j,bi,bj), st(i,j,bi,bj), ft(i,j,bi,bj),
U pHlocal(i,j), pCO2local(i,j), carbonate(i,j),
I i,j,k,bi,bj,myIter,myThid )
I i, j, k, bi, bj, myIter, myThid )
ENDDO
#ifdef CARBONCHEM_SOLVESAPHE
ENDIF
Expand Down Expand Up @@ -198,14 +198,16 @@ SUBROUTINE CALCITE_SATURATION( PTR_DIC, PTR_ALK, PTR_PO4,
C- else (case maskC = 0 , dry point)
omegaC(i,j,k,bi,bj) = 0. _d 0
ENDIF

C Update the pH3D array
pH3D(i,j,k,bi,bj) = pHlocal(i,j)
ENDDO
ENDDO

C Fill up some diagnostics of deep pH, pCO2, and carbonate local arrays
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(pHlocal ,'DIC3DPH ',k,1,2,bi,bj,myThid)
C Move outside i, j loops
C CALL DIAGNOSTICS_FILL(pHlocal ,'DIC3DPH ',k,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(pCO2local,'DIC3DPCO',k,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(carbonate,'DIC3DCO3',k,1,2,bi,bj,myThid)
ENDIF
Expand All @@ -216,6 +218,7 @@ SUBROUTINE CALCITE_SATURATION( PTR_DIC, PTR_ALK, PTR_PO4,

#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(pH3D ,'DIC3DPH ',0,Nr,1,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(omegaC ,'OMEGAC ',0,Nr,1,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(silicaDeep,'DIC3DSIT',0,Nr,1,bi,bj,myThid)
ENDIF
Expand Down
163 changes: 163 additions & 0 deletions pkg/dic/dic_calcite_read_pickup.F
Original file line number Diff line number Diff line change
@@ -0,0 +1,163 @@
#include "DIC_OPTIONS.h"

SUBROUTINE DIC_CALCITE_READ_PICKUP( pH3DisLoaded,
I myIter, myThid )

IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DIC_VARS.h"

C == Routine arguments ==
C myThid :: my Thread Id number
INTEGER myIter
INTEGER myThid

#ifdef ALLOW_DIC

C fn :: character buffer for creating filename
C prec :: precision of pickup files
C filePrec :: pickup-file precision (read from meta file)
C nbFields :: number of fields in pickup file (read from meta file)
C fldName :: Name of the field to read

C !FUNCTIONS
INTEGER MDS_RECLEN
EXTERNAL MDS_RECLEN

C !LOCAL VARIABLES:
C == Local variables ==
INTEGER i,j,k,bi,bj
INTEGER filePrec, ioUnit, prec, nbFields, nj
CHARACTER*(MAX_LEN_FNAM) fn, filNam
CHARACTER*(MAX_LEN_MBUF) msgBuf
CHARACTER*(8) fldName
CHARACTER*(8) suff
LOGICAL useCurrentDir, fileExist, pH3DisLoaded
INTEGER length_of_rec
_RL vec(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
CEOP
#ifdef DIC_CALCITE_SAT

pH3DisLoaded =.FALSE.

IF (pickupSuff.EQ.' ') THEN
IF ( rwSuffixType.EQ.0 ) THEN
WRITE(fn,'(A,I10.10)') 'pickup_dic_calcite.', myIter
ELSE
CALL RW_GET_SUFFIX( suff, startTime, myIter, myThid )
WRITE(fn,'(A,A)') 'pickup_dic_calcite.', suff
ENDIF
ELSE
WRITE(fn,'(A,A10)') 'pickup_dic_calcite.', pickupSuff
ENDIF
filePrec = precFloat64

C-- First check if pickup file exist
#ifdef ALLOW_MDSIO
useCurrentDir = .FALSE.
CALL MDS_CHECK4FILE(
I fn, '.data', 'DIC_CALCITE_READ_PICKUP',
O filNam, fileExist,
I useCurrentDir, myThid )
#else
STOP 'ABNORMAL END: S/R DIC_CALCITE_READ_PICKUP: Needs MDSIO pkg'
#endif
IF ( fileExist ) THEN
C-- Read pickup file

CALL READ_MFLDS_SET(
I fn,
O nbFields, filePrec,
I Nr, myIter, myThid )

_BEGIN_MASTER(myThid)
IF ( nbFields.GE.0 .AND. filePrec.NE.prec ) THEN
WRITE(msgBuf,'(2A,I4)') 'DIC_CALCITE_READ_PICKUP:',
& 'pickup-file binary precision do not match !'
CALL PRINT_ERROR( msgBuf, myThid )
WRITE(msgBuf,'(A,2(A,I4))') 'DIC_CALCITE_READ_PICKUP:',
& 'file prec.=', filePrec, ' but expecting prec.=', prec
CALL PRINT_ERROR( msgBuf, myThid )
STOP 'ABNORMAL END: S/R DIC_CALCITE_READ_PICKUP'
ENDIF
_END_MASTER( myThid )

IF ( nbFields.LE.0 ) THEN
C- No meta-file or old meta-file without List of Fields
ioUnit = errorMessageUnit
IF ( pickupStrictlyMatch ) THEN
WRITE(msgBuf,'(4A)') 'DIC_CALCITE_READ_PICKUP:',
& 'no field-list found in meta-file',
& ' => cannot check for strick-matching'
CALL PRINT_ERROR( msgBuf, myThid )
WRITE(msgBuf,'(4A)') 'DIC_CALCITE_READ_PICKUP:',
& 'try with " pickupStrictlyMatch=.FALSE.,"',
& ' in file: "data", NameList: "PARM03"'
CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
STOP 'ABNORMAL END: S/R DIC_CALCITE_READ_PICKUP'
ELSE
WRITE(msgBuf,'(4A)') 'WARNING >> DIC_CALCITE_READ_PICKUP:',
& ' no field-list found'
CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
IF ( nbFields.EQ.-1 ) THEN
C- No meta-file
WRITE(msgBuf,'(4A)') 'WARNING >> ',
& ' try to read pickup as currently written'
CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
ELSE
C- Old meta-file without List of Fields
WRITE(msgBuf,'(4A)') 'WARNING >> ',
& ' try to read pickup as it used to be written'
CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
WRITE(msgBuf,'(4A)') 'WARNING >> ',
& ' until checkpoint59l (2007 Dec 17)'
CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
ENDIF
ENDIF
ENDIF

nj = 0
C--- read pH 3-D fields for restart
fldName = 'pH3D '
CALL READ_MFLDS_3D_RL( fldName,
O vec,
& nj, prec, Nr, myIter, myThid )
CALL EXCH_3D_RL( vec, Nr, myThid )

DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO k=1,Nr
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
pH3D(i,j,k,bi,bj) = vec(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO

pH3DisLoaded =.TRUE.
ELSE
pH3DisLoaded = .FALSE.
ioUnit = errorMessageUnit

IF ( pickupStrictlyMatch ) THEN
WRITE(msgBuf,'(4A)') 'DIC_CALCITE_READ_PICKUP:',
& 'try with " pickupStrictlyMatch=.FALSE.,"',
& ' in file: "data", NameList: "PARM03"'
CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
STOP 'ABNORMAL END: S/R DIC_CALCITE_READ_PICKUP'
ELSE
WRITE(msgBuf,'(2A)') 'WARNING >> DIC_CALCITE_READ_PICKUP:',
& 'will restart from approximated pH'
CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
ENDIF
ENDIF
#endif /* DIC_CALCITE_SAT */
#endif /* ALLOW_DIC */

RETURN
END

0 comments on commit a7cfa71

Please sign in to comment.