Clean and fix pkg/dic initialisation (MITgcm#757)

* move 2 pkg/dic check & stop here

Better to have these in one place, and preferably checked early on at
the init-fixed stage, so move them from various places in pkg/dic src code.
In the future, more check & stop could be added here it needed.

* move pkg/dic check & stop to dic_init_fixed.F

* Clean and fix pkg/dic inititialisation

1) fix loading of forcing file for the case DIC_forcingCycle = 0
2) fix wrong/inconsistent timing params used for loading initial surface-silica file
3) move all pkg/dic initialise-varia level calls into dic_init_varia.F
   (previously most were called from gchem_init_vari.F), including
   DIC_READ_PICKUP and DIC_INI_ATMOS (previously called from dic_surfforcing_init.F)
4) move all forcing field initialisation to dic_ini_forcing.F
5) store logical var "pH_isLoaded" in common block ; this will allow to trigger
   a specific pH initialisation computation within the main pkg/dic calculation S/R
   that are called at each time-step
6) improve reading of "pickup_dic" file for future extension (e.g., adding an other field):
   read and process meta file with error/warning (depending on pickupStrictlyMatch = T/F)
   if missing file or fields.

* surrender to OpenAD

* Remove DIC_OPTIONS.h from top of DIC_VARS.h

and add it in the (only ?) file that does include DIC_VARS.h without DIC_OPTIONS.h

* adjusted to add back DIC_OPTIONS.h in module DIC_VARS

* Fix some comments and printed messages

* remove/move comments

* document improving pkg/dic initialisation

doc/tag-index

@@ -1,6 +1,13 @@
  Notes on tags used in MITgcmUV
  ==============================
 
+o pkg/dic:
+ - fix loading of forcing file for the case DIC_forcingCycle = 0 ;
+ - fix inconsistent timing params for loading initial surface-silica file ;
+ - clean-up initialise-varia level S/R, now all called from dic_init_varia.F
+ and move all forcing field initialisation to dic_ini_forcing.F ;
+ - read and process pickup_dic meta file with error/warning (depending on
+ pickupStrictlyMatch = T/F) if missing file or field.
 o pkg/ctrl:
  - reduce pkg/ctrl/optim.h to just "optimcycle" and rename into OPTIMCYCLE.h ;
  - keep optim_readparms.F namelist unchanged and define variables locally.

pkg/dic/DIC_VARS.h

@@ -1,6 +1,3 @@
-#ifdef ALLOW_OPENAD
-# include "DIC_OPTIONS.h"
-#endif
 C *==========================================================*
 C | DIC_VARS.h
 C | o Abiotic Carbon Variables
@@ -41,15 +38,17 @@ C =0 : Constant dissolution rate;
 C =1 : Follows (default) ;
 C =2 : Keir (1980) Geochem. Cosmochem. Acta. ;
 C =3 : Naviaux et al. 2019, Marine Chemistry
+C pH_isLoaded(1) :: = T when surface pH is loaded from pickup file
+C pH_isLoaded(2) :: = T when 3-D pH is loaded from pickup file
 
  COMMON /CARBON_NEEDS/
  & AtmospCO2, AtmosP, pH, pCO2, FluxCO2,
  & wind, fIce, Kwexch_Pre, silicaSurf,
  & calciteDissolRate, calciteDissolExp,
  & calcOmegaCalciteFreq, zca,
  & WsinkPIC, selectCalciteBottomRemin,
- & selectCalciteDissolution,
- & useCalciteSaturation, nIterCO3
+ & selectCalciteDissolution, nIterCO3,
+ & useCalciteSaturation, pH_isLoaded
 
  _RL AtmospCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
  _RL AtmosP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
@@ -68,6 +67,7 @@ C =3 : Naviaux et al. 2019, Marine Chemistry
  INTEGER selectCalciteDissolution
  INTEGER nIterCO3
  LOGICAL useCalciteSaturation
+ LOGICAL pH_isLoaded(2)
 
 #ifdef DIC_CALCITE_SAT
 C silicaDeep :: 3D-field of silicate concentration for pH and
@@ -413,4 +413,5 @@ C DIC_timeAve :: period over which DIC averages are calculated [s]
  _RL cfluxave (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
  _RL DIC_timeAve(nSx,nSy)
 #endif /* ALLOW_TIMEAVE */
+
 #endif /* DIC_BIOTIC */

pkg/dic/dic_fields_load.F

@@ -36,9 +36,6 @@ SUBROUTINE DIC_FIELDS_LOAD (
 #ifdef DIC_CALCITE_SAT
  INTEGER k
 #endif
-#ifdef READ_PAR
- CHARACTER*(MAX_LEN_MBUF) msgBuf
-#endif
 CEOP
 
  IF ( DIC_forcingCycle.GT.0. _d 0 ) THEN
@@ -258,12 +255,6 @@ SUBROUTINE DIC_FIELDS_LOAD (
  & + aWght*par1(i,j,bi,bj)
  ENDDO
  ENDDO
- ELSE
- WRITE(msgBuf,'(2A)')
- & ' DIC_FIELDS_LOAD: You need to provide ',
- & ' a file if you want to use READ_PAR'
- CALL PRINT_ERROR( msgBuf, myThid )
- STOP 'ABNORMAL END: S/R DIC_FIELDS_LOAD'
  ENDIF
 #endif
 #ifdef LIGHT_CHL

pkg/dic/dic_ini_forcing.F

@@ -26,43 +26,68 @@ SUBROUTINE DIC_INI_FORCING( myThid )
 
 #ifdef ALLOW_DIC
 
-c !LOCAL VARIABLES: ===================================================
- INTEGER bi,bj,i,j
+c !LOCAL VARIABLES: ====================================================
+ INTEGER bi, bj, i, j
+ INTEGER intimeP, intime0, intime1
+ _RL aWght, bWght
 #ifdef DIC_CALCITE_SAT
  INTEGER k
 #endif
-#if (defined (READ_PAR) && defined (USE_QSW))
- CHARACTER*(MAX_LEN_MBUF) msgBuf
-#endif
 
-C First call requires that we initialize everything to zero for safety
+C-- Initialise forcing arrays in common blocks
  DO bj = myByLo(myThid), myByHi(myThid)
  DO bi = myBxLo(myThid), myBxHi(myThid)
- DIC_ldRec(bi,bj) = 0
- ENDDO
- ENDDO
- CALL LEF_ZERO( dicwind0,myThid )
- CALL LEF_ZERO( dicwind1,myThid )
- CALL LEF_ZERO( atmosp0,myThid )
- CALL LEF_ZERO( atmosp1,myThid )
- CALL LEF_ZERO( silicaSurf0, myThid )
- CALL LEF_ZERO( silicaSurf1, myThid )
- CALL LEF_ZERO( ice0,myThid )
- CALL LEF_ZERO( ice1,myThid )
+
+C- Initialise DIC_VARS.h forcing-field arrays:
+ DO j=1-OLy,sNy+OLy
+ DO i=1-OLx,sNx+OLx
+ wind (i,j,bi,bj) = 0. _d 0
+ fIce (i,j,bi,bj) = 0. _d 0
+ AtmosP (i,j,bi,bj) = 0. _d 0
+ silicaSurf(i,j,bi,bj) = 0. _d 0
+#ifdef DIC_BIOTIC
+ par (i,j,bi,bj) = 0. _d 0
+ CHL (i,j,bi,bj) = 0. _d 0
+ InputFe (i,j,bi,bj) = 0. _d 0
+#endif
+ ENDDO
+ ENDDO
+#ifdef DIC_CALCITE_SAT
+ DO k=1,Nr
+ DO j=1-OLy,sNy+OLy
+ DO i=1-OLx,sNx+OLx
+ silicaDeep(i,j,k,bi,bj) = 0. _d 0
+ ENDDO
+ ENDDO
+ ENDDO
+#endif
+
+C- Initialise DIC_LOAD.h arrays:
+ DIC_ldRec(bi,bj) = 0
+ DO j=1-OLy,sNy+OLy
+ DO i=1-OLx,sNx+OLx
+ dicwind0(i,j,bi,bj) = 0. _d 0
+ dicwind1(i,j,bi,bj) = 0. _d 0
+ atmosp0 (i,j,bi,bj) = 0. _d 0
+ atmosp1 (i,j,bi,bj) = 0. _d 0
+ silicaSurf0(i,j,bi,bj) = 0. _d 0
+ silicaSurf1(i,j,bi,bj) = 0. _d 0
+ ice0 (i,j,bi,bj) = 0. _d 0
+ ice1 (i,j,bi,bj) = 0. _d 0
 #ifdef READ_PAR
- CALL LEF_ZERO( par0,myThid )
- CALL LEF_ZERO( par1,myThid )
+  par0 (i,j,bi,bj) = 0. _d 0
+  par1 (i,j,bi,bj) = 0. _d 0
 #endif
 #ifdef ALLOW_FE
- CALL LEF_ZERO( feinput0,myThid )
- CALL LEF_ZERO( feinput1,myThid )
+  feinput0(i,j,bi,bj) = 0. _d 0
+  feinput1(i,j,bi,bj) = 0. _d 0
 #endif
 #ifdef LIGHT_CHL
- CALL LEF_ZERO( chlinput,myThid )
+  chlinput(i,j,bi,bj) = 0. _d 0
 #endif
+ ENDDO
+ ENDDO
 #ifdef DIC_CALCITE_SAT
- 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
@@ -71,32 +96,24 @@ SUBROUTINE DIC_INI_FORCING( myThid )
  ENDDO
  ENDDO
  ENDDO
+#endif
+C- end bi,bj loops:
  ENDDO
  ENDDO
-#endif
 
-#ifdef READ_PAR
-#ifdef USE_QSW
- WRITE(msgBuf,'(2A)')
- & ' DIC_INI_FORCING: You can not use READ_PAR ',
- & ' and USE_QSW'
- CALL PRINT_ERROR( msgBuf, myThid )
- STOP 'ABNORMAL END: S/R DIC_INI_FORCING'
-#endif
-#endif
+C ======================================================================
 
-C set reasonable values to those that need at least something
-  DO bj = myByLo(myThid), myByHi(myThid)
-  DO bi = myBxLo(myThid), myBxHi(myThid)
+C-- Set reasonable values to those that need at least something
+ DO bj = myByLo(myThid), myByHi(myThid)
+ DO bi = myBxLo(myThid), myBxHi(myThid)
  DO j=1-OLy,sNy+OLy
  DO i=1-OLx,sNx+OLx
- WIND(i,j,bi,bj) = 5. _d 0*maskC(i,j,1,bi,bj)
+ wind(i,j,bi,bj) = 5. _d 0*maskC(i,j,1,bi,bj)
  AtmosP(i,j,bi,bj) = 1. _d 0*maskC(i,j,1,bi,bj)
  silicaSurf(i,j,bi,bj) = 7.6838 _d -3*maskC(i,j,1,bi,bj)
  fIce(i,j,bi,bj) = 0. _d 0
- FluxCO2(i,j,bi,bj) = 0. _d 0
 #ifdef READ_PAR
- PAR(i,j,bi,bj) = 100. _d 0*maskC(i,j,1,bi,bj)
+ par(i,j,bi,bj) = 100. _d 0*maskC(i,j,1,bi,bj)
 #endif
 #ifdef LIGHT_CHL
 C If the chlorophyll climatology is not provided, use this default value.
@@ -119,8 +136,131 @@ SUBROUTINE DIC_INI_FORCING( myThid )
  ENDIF
 #endif
 C- end bi,bj loops
- ENDDO
  ENDDO
+ ENDDO
+
+C ======================================================================
+
+C-- Load constant forcing-file (used if DIC_forcingCycle=0 )
+c IF ( DIC_forcingCycle .LE. zeroRL ) THEN
+ IF ( DIC_windFile .NE. ' ' ) THEN
+ CALL READ_REC_XY_RL( DIC_windFile, wind,
+ & 1, nIter0, myThid )
+ _EXCH_XY_RL( wind, myThid )
+ ENDIF
+ IF ( DIC_atmospFile .NE. ' ' ) THEN
+ CALL READ_REC_XY_RL( DIC_atmospFile, AtmosP,
+ & 1, nIter0, myThid )
+ _EXCH_XY_RL( AtmosP, myThid )
+ ENDIF
+ IF ( DIC_silicaFile .NE. ' ' ) THEN
+ CALL READ_REC_XY_RL( DIC_silicaFile, silicaSurf,
+ & 1, nIter0, myThid )
+ _EXCH_XY_RL( silicaSurf, myThid )
+ ENDIF
+ IF ( DIC_iceFile .NE. ' ' ) THEN
+ CALL READ_REC_XY_RL( DIC_iceFile, fIce,
+ & 1, nIter0, myThid )
+ _EXCH_XY_RL( fIce, myThid )
+ ENDIF
+#ifdef DIC_BIOTIC
+c#ifdef READ_PAR
+ IF ( DIC_parFile .NE. ' ' ) THEN
+ CALL READ_REC_XY_RL( DIC_parFile, par,
+ & 1, nIter0, myThid )
+ _EXCH_XY_RL( par, myThid )
+ ENDIF
+c#endif
+c#ifdef LIGHT_CHL
+C-- Load chlorophyll climatology data, unit for chlorophyll : mg/m3
+ IF ( DIC_chlaFile .NE. ' ' ) THEN
+ CALL READ_REC_XY_RL( DIC_chlaFile, CHL,
+ & 1, nIter0, myThid )
+ _EXCH_XY_RL( CHL, myThid )
+ ENDIF
+c#endif
+c#ifdef ALLOW_FE
+ IF ( DIC_ironFile .NE. ' ' ) THEN
+ CALL READ_REC_XY_RL( DIC_ironFile, InputFe,
+ & 1, nIter0, myThid )
+ _EXCH_XY_RL( InputFe, myThid )
+ ENDIF
+c#endif
+#endif /* DIC_BIOTIC */
+#ifdef DIC_CALCITE_SAT
+ IF ( useCalciteSaturation .AND.
+ & DIC_deepSilicaFile .NE. ' ' ) THEN
+ CALL READ_REC_XYZ_RL( DIC_deepSilicaFile, silicaDeep,
+ & 1, nIter0, myThid )
+ _EXCH_XYZ_RL( silicaDeep, myThid )
+ ENDIF
+#endif
+c ENDIF
+
+C ======================================================================
+
+ IF ( DIC_forcingCycle .GT. zeroRL ) THEN
+C-- Load some initial forcing from file
+
+C- Read in surface silica field (used to compute initial surf. pH)
+C Note: For this purpose, might be accurate enough to use first record
+C in silica file (already loaded) and skip all this block.
+
+C- Get reccord number and weight for time interpolation:
+ CALL GET_PERIODIC_INTERVAL(
+ O intimeP, intime0, intime1, bWght, aWght,
+ I DIC_forcingCycle, DIC_forcingPeriod,
+ I deltaTClock, startTime, myThid )
+
+ _BARRIER
+ _BEGIN_MASTER(myThid)
+ WRITE(standardMessageUnit,'(A,I10,A,2(2I5,A))')
+ & ' DIC_INI_FORCING, it=', nIter0,
+ & ' : Reading new data, i0,i1=', intime0, intime1
+ _END_MASTER(myThid)
+
+ IF ( DIC_silicaFile .NE. ' ' ) THEN
+C- If file provided for surface silicate, read it in
+ CALL READ_REC_XY_RS( DIC_silicaFile,silicaSurf0,intime0,
+ & nIter0,myThid )
+ CALL READ_REC_XY_RS( DIC_silicaFile,silicaSurf1,intime1,
+ & nIter0,myThid )
+ DO bj = myByLo(myThid), myByHi(myThid)
+ DO bi = myBxLo(myThid), myBxHi(myThid)
+ DO j=1-OLy,sNy+OLy
+ DO i=1-OLx,sNx+OLx
+ silicaSurf(i,j,bi,bj) = bWght*silicaSurf0(i,j,bi,bj)
+ & + aWght*silicaSurf1(i,j,bi,bj)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+#ifndef ALLOW_AUTODIFF
+#ifdef DIC_CALCITE_SAT
+ ELSEIF ( DIC_deepSilicaFile .NE. ' ' ) THEN
+C- If no surface silicate file but deep (3d) silicate provided, use top level
+ k = 1
+ CALL READ_REC_XYZ_RS( DIC_deepSilicaFile, silicaDeep0,
+ & intime0, nIter0, myThid )
+ CALL READ_REC_XYZ_RS( DIC_deepSilicaFile, silicaDeep1,
+ & intime1, nIter0, myThid )
+ DO bj = myByLo(myThid), myByHi(myThid)
+ DO bi = myBxLo(myThid), myBxHi(myThid)
+ DO j=1-OLy,sNy+OLy
+ DO i=1-OLx,sNx+OLx
+ silicaSurf(i,j,bi,bj) = bWght*silicaDeep0(i,j,k,bi,bj)
+ & + aWght*silicaDeep1(i,j,k,bi,bj)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+#endif /* DIC_CALCITE_SAT */
+#endif /* ndef ALLOW_AUTODIFF */
+ ENDIF
+ _EXCH_XY_RL( silicaSurf, myThid )
+
+C end if DIC_forcingCycle > 0
+ ENDIF
 
 #endif /* ALLOW_DIC */
  RETURN