From 51e381e9c9a1bef17c6feb83a95af189631764de Mon Sep 17 00:00:00 2001 From: Jean-Michel Campin Date: Tue, 7 Nov 2023 13:00:07 -0500 Subject: [PATCH] Clean and fix pkg/dic initialisation (#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 | 7 + pkg/dic/DIC_VARS.h | 11 +- pkg/dic/dic_fields_load.F | 9 - pkg/dic/dic_ini_forcing.F | 222 +++++++++++++++++++----- pkg/dic/dic_init_fixed.F | 32 +++- pkg/dic/dic_init_varia.F | 56 ++++-- pkg/dic/dic_read_pickup.F | 249 ++++++++++++++++++++++----- pkg/dic/dic_surfforcing_init.F | 118 +------------ pkg/gchem/gchem_init_vari.F | 6 - pkg/openad/inner_do_loop.F | 5 +- tools/OAD_support/cb2mGetModules.csh | 2 +- 11 files changed, 487 insertions(+), 230 deletions(-) diff --git a/doc/tag-index b/doc/tag-index index fac66c9415..21fc5c3a49 100644 --- a/doc/tag-index +++ b/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. diff --git a/pkg/dic/DIC_VARS.h b/pkg/dic/DIC_VARS.h index 0f768b4a66..6028857397 100644 --- a/pkg/dic/DIC_VARS.h +++ b/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,6 +38,8 @@ 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, @@ -48,8 +47,8 @@ C =3 : Naviaux et al. 2019, Marine Chemistry & 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 */ diff --git a/pkg/dic/dic_fields_load.F b/pkg/dic/dic_fields_load.F index e154c9b466..6c41510bea 100644 --- a/pkg/dic/dic_fields_load.F +++ b/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 diff --git a/pkg/dic/dic_ini_forcing.F b/pkg/dic/dic_ini_forcing.F index fd3be33b7e..28daf93a41 100644 --- a/pkg/dic/dic_ini_forcing.F +++ b/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 diff --git a/pkg/dic/dic_init_fixed.F b/pkg/dic/dic_init_fixed.F index c494f90a3f..ede4305742 100644 --- a/pkg/dic/dic_init_fixed.F +++ b/pkg/dic/dic_init_fixed.F @@ -27,7 +27,7 @@ SUBROUTINE DIC_INIT_FIXED( myThid ) #ifdef ALLOW_DIC INTEGER k INTEGER iUnit -#ifdef DIC_BIOTIC +#if ( defined DIC_BIOTIC || defined READ_PAR ) CHARACTER*(MAX_LEN_MBUF) msgBuf #endif @@ -35,7 +35,7 @@ SUBROUTINE DIC_INIT_FIXED( myThid ) _BEGIN_MASTER(myThid) -C set up coefficients for DIC chemistry +C-- Set up coefficients for DIC chemistry C define Schmidt no. coefficients for CO2 sca1 = 2073.1 _d 0 sca2 = -125.62 _d 0 @@ -61,7 +61,7 @@ SUBROUTINE DIC_INIT_FIXED( myThid ) oB3= -8.17083 _d -3 oC0= -4.88682 _d -7 -C Set other constant/flag +C-- Set other constant/flag IF ( dic_int1.EQ.2 ) THEN CALL MDSFINDUNIT( iUnit, myThid ) @@ -93,12 +93,36 @@ SUBROUTINE DIC_INIT_FIXED( myThid ) WRITE(msgBuf,'(A)') '// DIC_INIT_FIXED parameters :' CALL PRINT_MESSAGE(msgBuf,iUnit,SQUEEZE_RIGHT,myThid) - CALL WRITE_0D_I( nlev, INDEX_NONE, 'nlev =', + CALL WRITE_0D_I( nlev, INDEX_NONE, 'nlev =', & ' /* Number of level over which Bio act is computed */') #endif /* DIC_BIOTIC */ _END_MASTER(myThid) +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| +C-- Check some option & parameter combinations: + +C for now, just collect here few check & stop from various places: +#ifdef READ_PAR +#ifdef USE_QSW +c WRITE(msgBuf,'(2A)') ' DIC_INI_FORCING: ', + WRITE(msgBuf,'(2A)') 'DIC_INIT_FIXED: ', + & 'You can not use READ_PAR and USE_QSW together' + CALL PRINT_ERROR( msgBuf, myThid ) + STOP 'ABNORMAL END: S/R DIC_INIT_FIXED' +#endif + IF ( DIC_forcingCycle.GT.0. _d 0 + & .AND. DIC_parFile .EQ. ' ' ) THEN +c WRITE(msgBuf,'(2A)') ' DIC_FIELDS_LOAD: ', + WRITE(msgBuf,'(2A)') 'DIC_INIT_FIXED: ', + & 'You need to provide a file if you want to use READ_PAR' + CALL PRINT_ERROR( msgBuf, myThid ) + STOP 'ABNORMAL END: S/R DIC_INIT_FIXED' + ENDIF +#endif /* READ_PAR */ + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| + #ifdef ALLOW_MNC IF ( useMNC ) THEN CALL DIC_MNC_INIT( myThid ) diff --git a/pkg/dic/dic_init_varia.F b/pkg/dic/dic_init_varia.F index 1e9d9b2fd9..383832951c 100644 --- a/pkg/dic/dic_init_varia.F +++ b/pkg/dic/dic_init_varia.F @@ -1,4 +1,5 @@ #include "DIC_OPTIONS.h" +#undef USE_OLD_READ_PICKUP C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP @@ -18,27 +19,29 @@ SUBROUTINE DIC_INIT_VARIA( myThid ) #include "GRID.h" #include "DIC_VARS.h" #include "DIC_ATMOS.h" +#include "PTRACERS_SIZE.h" +#include "PTRACERS_PARAMS.h" #ifdef ALLOW_COST # include "DIC_COST.h" #endif C !INPUT PARAMETERS: -C myThid :: thread number +C myThid :: my Thread Id number INTEGER myThid CEOP #ifdef ALLOW_DIC -#ifdef DIC_BIOTIC INTEGER i,j, bi,bj # ifdef DIC_CALCITE_SAT INTEGER k # endif -#endif c CHARACTER*(MAX_LEN_MBUF) msgBuf C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Initialise variable in common block DIC_ATMOS +#ifndef ALLOW_AUTODIFF _BEGIN_MASTER(myThid) +#endif total_atmos_carbon = 0. _d 0 total_ocean_carbon = 0. _d 0 total_atmos_carbon_year = 0. _d 0 @@ -49,31 +52,62 @@ SUBROUTINE DIC_INIT_VARIA( myThid ) #ifdef ALLOW_COST totcost = 0. _d 0 #endif + pH_isLoaded(1) = .FALSE. + pH_isLoaded(2) = .FALSE. +#ifndef ALLOW_AUTODIFF _END_MASTER(myThid) _BARRIER +#endif -#ifdef DIC_BIOTIC -C-- Initialise alpha & rain_ratio fields with fixed (& Uniform) values +C-- Initialise variables other than forcing (done in ini_forcing): DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy - DO i=1-OLx,sNx+OLx - alpha(i,j,bi,bj) = alphaUniform - rain_ratio(i,j,bi,bj) = rainRatioUniform - ENDDO + DO i=1-OLx,sNx+OLx + FluxCO2(i,j,bi,bj) = 0. _d 0 + pH(i,j,bi,bj) = 8. _d 0 + ENDDO + ENDDO +#ifdef DIC_BIOTIC +C- Initialise alpha & rain_ratio fields with fixed (& Uniform) values + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + alpha(i,j,bi,bj) = alphaUniform + rain_ratio(i,j,bi,bj) = rainRatioUniform + ENDDO ENDDO +#endif /* DIC_BIOTIC */ #ifdef DIC_CALCITE_SAT DO k = 1, Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx - omegaC(i,j,k,bi,bj) = 0. _d 0 + omegaC(i,j,k,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO #endif /* DIC_CALCITE_SAT */ ENDDO ENDDO -#endif /* DIC_BIOTIC */ + +C-- Initialise forcing field variables (was called from gchem_init_vari.F) + CALL DIC_INI_FORCING(myThid) + +C- read pickup + IF ( nIter0.GT.PTRACERS_Iter0 .OR. + & (nIter0.EQ.PTRACERS_Iter0 .AND. pickupSuff.NE.' ') + & ) THEN +C Read pH from a pickup file if needed + CALL DIC_READ_PICKUP( nIter0, myThid ) + ENDIF + +C- Move these S/R call here (were previously in gchem_init_vari.F), except +C DIC_INI_ATMOS which was called from the top of DIC_SURFFORCING_INIT + CALL DIC_INI_ATMOS( startTime, nIter0, myThid ) + CALL DIC_SURFFORCING_INIT(myThid) + CALL DIC_BIOTIC_INIT(myThid) +# ifdef ALLOW_CTRL + CALL DIC_SET_CONTROL(myThid) +# endif #endif /* ALLOW_DIC */ diff --git a/pkg/dic/dic_read_pickup.F b/pkg/dic/dic_read_pickup.F index 9fbd4e7b36..9a2eb235c8 100644 --- a/pkg/dic/dic_read_pickup.F +++ b/pkg/dic/dic_read_pickup.F @@ -1,85 +1,252 @@ #include "DIC_OPTIONS.h" - SUBROUTINE DIC_READ_PICKUP( - O pH_isLoaded, - I myIter, myThid ) +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| +CBOP +C !ROUTINE: DIC_READ_PICKUP +C !INTERFACE: + SUBROUTINE DIC_READ_PICKUP( myIter, myThid ) + +C !DESCRIPTION: +C Reads current state of DIC from a pickup file + +C !USES: IMPLICIT NONE -C === Global variables === +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 - LOGICAL pH_isLoaded +C !INPUT PARAMETERS: +C myIter :: time-step number in simulation +C myThid :: my Thread Id number INTEGER myIter INTEGER myThid #ifdef ALLOW_DIC #ifdef DIC_BIOTIC - -C !FUNCTIONS +C !FUNCTIONS: + INTEGER ILNBLNK + EXTERNAL ILNBLNK C !LOCAL VARIABLES: +C fn :: character buffer for creating filename +C fp :: 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 missFldList :: List of missing fields (attempted to read but not found) +C missFldDim :: Dimension of missing fields list array: missFldList +C nMissing :: Number of missing fields (attempted to read but not found) +C j :: loop index +C nj :: record number +C ioUnit :: temp for writing msg unit +C msgBuf :: Informational/error message buffer + INTEGER fp + INTEGER filePrec, nbFields + INTEGER missFldDim, nMissing + INTEGER j, nj, ioUnit + PARAMETER( missFldDim = 12 ) CHARACTER*(10) suff - CHARACTER*(MAX_LEN_FNAM) fn, filNam + CHARACTER*(MAX_LEN_FNAM) fn + CHARACTER*(8) missFldList(missFldDim) CHARACTER*(MAX_LEN_MBUF) msgBuf - LOGICAL useCurrentDir, fileExist - INTEGER fp, ioUnit + CHARACTER*(MAX_LEN_FNAM) tmpNam +C- note: to avoid beeing caught by tools/OAD_support/stop2print.sed, +C change "stopFlag" to "StopFlag" + LOGICAL useCurrentDir, fileExist, StopFlag + INTEGER iL CEOP - pH_isLoaded =.FALSE. - ioUnit = errorMessageUnit +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| -C-- pickup file name : - IF (pickupSuff.EQ.' ') THEN + IF ( pickupSuff.EQ.' ' ) THEN IF ( rwSuffixType.EQ.0 ) THEN WRITE(fn,'(A,I10.10)') 'pickup_dic.', myIter ELSE CALL RW_GET_SUFFIX( suff, startTime, myIter, myThid ) WRITE(fn,'(A,A)') 'pickup_dic.', suff ENDIF - ELSE + ELSE WRITE(fn,'(A,A10)') 'pickup_dic.', pickupSuff - ENDIF - fp = precFloat64 + ENDIF + fp = precFloat64 -C-- First check if pickup file exist + CALL READ_MFLDS_SET( + I fn, + O nbFields, filePrec, + I Nr, myIter, myThid ) + _BEGIN_MASTER( myThid ) +c IF ( filePrec.NE.0 .AND. filePrec.NE.fp ) THEN + IF ( nbFields.GE.0 .AND. filePrec.NE.fp ) THEN + WRITE(msgBuf,'(2A,I4)') 'DIC_READ_PICKUP: ', + & 'pickup-file binary precision do not match !' + CALL PRINT_ERROR( msgBuf, myThid ) + WRITE(msgBuf,'(A,2(A,I4))') 'DIC_READ_PICKUP: ', + & 'file prec.=', filePrec, ' but expecting prec.=', fp + CALL PRINT_ERROR( msgBuf, myThid ) + CALL ALL_PROC_DIE( 0 ) + STOP 'ABNORMAL END: S/R DIC_READ_PICKUP (data-prec Pb)' + ENDIF + _END_MASTER( myThid ) + + ioUnit = errorMessageUnit + StopFlag = .FALSE. + IF ( nbFields.LE.0 ) THEN +C- No meta-file or old meta-file without List of Fields + IF ( pickupStrictlyMatch ) THEN + WRITE(msgBuf,'(4A)') 'DIC_READ_PICKUP: ', + & 'no field-list found in meta-file', + & ' => cannot check for strick-matching' + CALL PRINT_ERROR( msgBuf, myThid ) + WRITE(msgBuf,'(4A)') 'DIC_READ_PICKUP: ', + & 'try with " pickupStrictlyMatch=.FALSE.,"', + & ' in file: "data", NameList: "PARM03"' + CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid ) + StopFlag = .TRUE. + ELSE + WRITE(msgBuf,'(4A)') 'WARNING >> DIC_READ_PICKUP: ', + & ' no field-list found' + CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid ) + IF ( nbFields.EQ.-1 ) THEN +C- No meta-file: then check if binary pickup file (i.e., ".data") exist #ifdef ALLOW_MDSIO - useCurrentDir = .FALSE. - CALL MDS_CHECK4FILE( + useCurrentDir = .FALSE. + CALL MDS_CHECK4FILE( I fn, '.data', 'DIC_READ_PICKUP', - O filNam, fileExist, + O tmpNam, fileExist, I useCurrentDir, myThid ) #else - STOP 'ABNORMAL END: S/R DIC_READ_PICKUP: Needs MDSIO pkg' + STOP 'ABNORMAL END: S/R DIC_READ_PICKUP: Needs MDSIO pkg' #endif + IF ( fileExist ) THEN + WRITE(msgBuf,'(4A)') 'WARNING >> ', + & ' try to read pickup as currently written' + CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid ) + ELSE + iL = ILNBLNK(fn) + WRITE(msgBuf,'(4A)') 'DIC_READ_PICKUP: ', + & 'missing both "meta" & "data" files for "', fn(1:iL), '"' + CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid ) + nbFields = -2 + ENDIF + ELSE +C- Old meta-file without List of Fields +c WRITE(msgBuf,'(4A)') 'WARNING >> ', +c & ' try to read pickup as it used to be written' +c CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid ) +c WRITE(msgBuf,'(4A)') 'WARNING >> ', +c & ' until checkpoint59l (2007 Dec 17)' +c CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid ) + WRITE(msgBuf,'(4A)') 'DIC_READ_PICKUP: ', + & 'no field-list found in meta-file' + CALL PRINT_ERROR( msgBuf, myThid ) + StopFlag = .TRUE. + ENDIF + ENDIF + ENDIF + IF ( StopFlag ) THEN + CALL ALL_PROC_DIE( myThid ) + STOP 'ABNORMAL END: S/R DIC_READ_PICKUP' + ENDIF - IF ( fileExist ) THEN -C-- Read pickup file - CALL READ_REC_3D_RL( fn, fp, 1, pH, 1, myIter, myThid ) - pH_isLoaded = .TRUE. +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - _EXCH_XY_RL( pH, myThid ) - ELSE - pH_isLoaded = .FALSE. - IF ( pickupStrictlyMatch ) THEN + IF ( nbFields.EQ.0 ) THEN +C--- Old meta-file without List of Fields: use the old way to read pickup + + ELSEIF ( nbFields.NE.-2 ) THEN +C--- New way to read DIC pickup: + nj = 0 +C--- read DIC 3-D fields for restart +#ifdef DIC_CALCITE_SAT + IF ( useCalciteSaturation ) THEN +c CALL READ_MFLDS_3D_RL( 'DIC_pH3d', pH3D, +c & nj, fp, Nr, myIter, myThid ) +c _BEGIN_MASTER( myThid ) +c pH_isLoaded(2) = .TRUE. +c _END_MASTER( myThid ) + ENDIF +#endif + +C- switch to 2-D fields: + nj = nj*Nr +C--- read DIC 2-D fields for restart + CALL READ_MFLDS_3D_RL( 'DIC_pH2d', pH, + & nj, fp, 1 , myIter, myThid ) + _BEGIN_MASTER( myThid ) + pH_isLoaded(1) = .TRUE. + _END_MASTER( myThid ) + +C-- end: new way to read pickup file + ENDIF + +C-- Check for missing fields: + nMissing = missFldDim + CALL READ_MFLDS_CHECK( + O missFldList, + U nMissing, + I myIter, myThid ) + IF ( nMissing.GT.missFldDim ) THEN + WRITE(msgBuf,'(2A,I4)') 'DIC_READ_PICKUP: ', + & 'missing fields list has been truncated to', missFldDim + CALL PRINT_ERROR( msgBuf, myThid ) + CALL ALL_PROC_DIE( myThid ) + STOP 'ABNORMAL END: S/R DIC_READ_PICKUP (list-size Pb)' + ENDIF + IF ( nMissing.GE.1 ) THEN + DO j=1,nMissing + IF ( missFldList(nj) .EQ. 'DIC_pH2d' ) THEN + _BEGIN_MASTER( myThid ) + pH_isLoaded(1) = .FALSE. + _END_MASTER( myThid ) + ELSEIF ( missFldList(nj) .EQ. 'DIC_pH3d' ) THEN + _BEGIN_MASTER( myThid ) + pH_isLoaded(2) = .FALSE. + _END_MASTER( myThid ) + ELSE + StopFlag = .TRUE. + WRITE(msgBuf,'(4A)') 'DIC_READ_PICKUP: ', + & 'cannot restart without field "',missFldList(nj),'"' + CALL PRINT_ERROR( msgBuf, myThid ) + ENDIF + ENDDO + IF ( pickupStrictlyMatch .AND. .NOT.StopFlag ) THEN + StopFlag = .TRUE. WRITE(msgBuf,'(4A)') 'DIC_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_READ_PICKUP' - ELSE - WRITE(msgBuf,'(2A)') 'WARNING >> DIC_READ_PICKUP: ', - & 'will restart from approximated pH' - CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid ) + CALL PRINT_ERROR( msgBuf, myThid ) ENDIF - ENDIF + ENDIF + IF ( StopFlag ) THEN + CALL ALL_PROC_DIE( myThid ) + STOP 'ABNORMAL END: S/R DIC_READ_PICKUP' + ENDIF + + _BEGIN_MASTER( myThid ) + IF ( .NOT.pH_isLoaded(1) ) THEN + WRITE(msgBuf,'(2A)') 'WARNING >> DIC_READ_PICKUP: ', + & 'will restart from approximated 2-D pH' + CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid ) + ENDIF + IF ( useCalciteSaturation .AND. .NOT.pH_isLoaded(2) ) THEN + WRITE(msgBuf,'(2A)') 'WARNING >> DIC_READ_PICKUP: ', + & 'will restart from approximated 3-D pH' +c CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid ) + ENDIF + _END_MASTER( myThid ) + +C-- Update overlap regions: + CALL EXCH_XY_RL( pH, myThid ) +#ifdef DIC_CALCITE_SAT + IF ( useCalciteSaturation ) THEN +c CALL EXCH_3D_RL( pH3D, Nr, myThid ) + ENDIF +#endif -#endif /* DIC_BIOTIC */ -#endif /* ALLOW_DIC */ +#endif /* DIC_BIOTIC */ +#endif /* ALLOW_DIC */ RETURN END diff --git a/pkg/dic/dic_surfforcing_init.F b/pkg/dic/dic_surfforcing_init.F index 14e944d95c..a359cf6c90 100644 --- a/pkg/dic/dic_surfforcing_init.F +++ b/pkg/dic/dic_surfforcing_init.F @@ -23,7 +23,6 @@ SUBROUTINE DIC_SURFFORCING_INIT( #include "PTRACERS_SIZE.h" #include "PTRACERS_PARAMS.h" #include "PTRACERS_FIELDS.h" -#include "DIC_LOAD.h" C !INPUT PARAMETERS: =================================================== C myThid :: thread number @@ -32,13 +31,10 @@ SUBROUTINE DIC_SURFFORCING_INIT( #ifdef ALLOW_DIC C !LOCAL VARIABLES: ==================================================== - INTEGER i,j, kLev, it - INTEGER intimeP, intime0, intime1 - _RL aWght, bWght, co3dummy -C Number of iterations for pCO2 solvers... -C Solubility relation coefficients -C local variables for carbon chem + INTEGER i,j, kLev, it INTEGER iMin,iMax,jMin,jMax, bi, bj + _RL co3dummy +C local variables for carbon chem _RL surfalk(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL surfphos(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL surfsi(1-OLx:sNx+OLx,1-OLy:sNy+OLy) @@ -47,109 +43,21 @@ SUBROUTINE DIC_SURFFORCING_INIT( _RL surfdic(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER iprt,jprt LOGICAL debugPrt - LOGICAL pH_isLoaded +#ifdef ALLOW_DEBUG C msgBuf :: Informational/error message buffer CHARACTER*(MAX_LEN_MBUF) msgBuf +#endif CEOP - kLev = 1 - - CALL DIC_INI_ATMOS( startTime, nIter0, myThid ) - -C ================================================================= - - IF ( periodicExternalForcing ) THEN - -C read in surface silica field -C-- Get reccord number and weight for time interpolation: - CALL GET_PERIODIC_INTERVAL( - O intimeP, intime0, intime1, bWght, aWght, - I externForcingCycle, externForcingPeriod, - I deltaTClock, startTime, myThid ) - - _BARRIER - _BEGIN_MASTER(myThid) - WRITE(standardMessageUnit,'(A,I10,A,2(2I5,A))') - & ' DIC_SURFFORCING_INIT, it=', nIter0, - & ' : Reading new data, i0,i1=', intime0, intime1 - _END_MASTER(myThid) - - IF ( DIC_silicaFile .NE. ' ' ) THEN - CALL READ_REC_XY_RS( DIC_silicaFile,silicaSurf0,intime0, - & nIter0,myThid ) - CALL READ_REC_XY_RS( DIC_silicaFile,silicaSurf1,intime1, - & nIter0,myThid ) -#ifndef ALLOW_AUTODIFF -#ifdef DIC_CALCITE_SAT - ELSEIF ( DIC_deepSilicaFile .NE. ' ' ) THEN - CALL READ_REC_XYZ_RS( DIC_deepSilicaFile, silicaDeep0, - & intime0, nIter0, myThid ) - CALL READ_REC_XYZ_RS( DIC_deepSilicaFile, silicaDeep1, - & intime1, nIter0, myThid ) -#endif /* DIC_CALCITE_SAT */ -#endif /* ndef ALLOW_AUTODIFF */ - ENDIF - -c#ifdef ALLOW_OFFLINE -C-- This call has been moved to S/R OFFLINE_INIT_VARIA -c IF ( useOffLine ) THEN -c CALL OFFLINE_FIELDS_LOAD( startTime, nIter0, myThid ) -c ENDIF -c#endif - - _EXCH_XY_RS(silicaSurf0, myThid ) - _EXCH_XY_RS(silicaSurf1, myThid ) - - IF ( DIC_silicaFile .NE. ' ' ) THEN - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO j=1-OLy,sNy+OLy - DO i=1-OLx,sNx+OLx -C If file provided for surface silicate, read it in - 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 - DO bj = myByLo(myThid), myByHi(myThid) - DO bi = myBxLo(myThid), myBxHi(myThid) - DO j=1-OLy,sNy+OLy - DO i=1-OLx,sNx+OLx -C If no surface silicate file but deep (3d) silicate provided, use top level - silicaSurf(i,j,bi,bj) = bWght*silicaDeep0(i,j,kLev,bi,bj) - & + aWght*silicaDeep1(i,j,kLev,bi,bj) - ENDDO - ENDDO - ENDDO - ENDDO -#endif /* DIC_CALCITE_SAT */ -#endif /* ndef ALLOW_AUTODIFF */ - ENDIF - -C end periodicExternalForcing - ENDIF - C ================================================================= + kLev = 1 jMin=1 jMax=sNy iMin=1 iMax=sNx - DO bj=myByLo(myThid),myByHi(myThid) - DO bi=myBxLo(myThid),myBxHi(myThid) - DO j=1-OLy,sNy+OLy - DO i=1-OLx,sNx+OLx - pH(i,j,bi,bj) = 8. _d 0 - ENDDO - ENDDO - ENDDO - ENDDO - +C Solubility relation coefficients DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy @@ -174,16 +82,6 @@ SUBROUTINE DIC_SURFFORCING_INIT( ENDDO ENDDO - pH_isLoaded = .FALSE. - IF ( nIter0.GT.PTRACERS_Iter0 .OR. - & (nIter0.EQ.PTRACERS_Iter0 .AND. pickupSuff.NE.' ') - & ) THEN -C Read pH from a pickup file if needed - CALL DIC_READ_PICKUP( - O pH_isLoaded, - I nIter0, myThid ) - ENDIF - DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) @@ -246,7 +144,7 @@ SUBROUTINE DIC_SURFFORCING_INIT( C==================================================================== - IF ( .NOT.pH_isLoaded ) THEN + IF ( .NOT.pH_isLoaded(1) ) THEN C set guess of pH for first step here #ifdef ALLOW_DEBUG IF (debugMode) THEN diff --git a/pkg/gchem/gchem_init_vari.F b/pkg/gchem/gchem_init_vari.F index 38d53be326..8fd5c45d5f 100644 --- a/pkg/gchem/gchem_init_vari.F +++ b/pkg/gchem/gchem_init_vari.F @@ -67,12 +67,6 @@ SUBROUTINE GCHEM_INIT_VARI(myThid ) IF ( useDIC ) THEN # endif /* ALLOW_AUTODIFF */ CALL DIC_INIT_VARIA(myThid) - CALL DIC_INI_FORCING(myThid) - CALL DIC_SURFFORCING_INIT(myThid) - CALL DIC_BIOTIC_INIT(myThid) -# ifdef ALLOW_CTRL - CALL DIC_SET_CONTROL(myThid) -# endif # ifndef ALLOW_AUTODIFF ENDIF # endif diff --git a/pkg/openad/inner_do_loop.F b/pkg/openad/inner_do_loop.F index 50b6af2e3a..4d75488c5b 100644 --- a/pkg/openad/inner_do_loop.F +++ b/pkg/openad/inner_do_loop.F @@ -1,5 +1,3 @@ -#include "PACKAGES_CONFIG.h" -#include "CPP_OPTIONS.h" #include "OPENAD_OPTIONS.h" c#ifdef ALLOW_AUTODIFF c# include "AUTODIFF_OPTIONS.h" @@ -13,6 +11,9 @@ #ifdef ALLOW_GENERIC_ADVDIFF # include "GAD_OPTIONS.h" #endif +#ifdef ALLOW_DIC +# include "DIC_OPTIONS.h" +#endif #ifdef ALLOW_GMREDI # include "GMREDI_OPTIONS.h" #endif diff --git a/tools/OAD_support/cb2mGetModules.csh b/tools/OAD_support/cb2mGetModules.csh index 5b44a9c8ee..ab1173f149 100755 --- a/tools/OAD_support/cb2mGetModules.csh +++ b/tools/OAD_support/cb2mGetModules.csh @@ -16,7 +16,7 @@ echo '#endif' >> ${fileName}_temp echo '#ifdef ALLOW_ECCO' >> ${fileName}_temp echo '# include "ECCO_OPTIONS.h"' >> ${fileName}_temp echo '#endif' >> ${fileName}_temp -if ( ${fileName} == 'DIC_LOAD' ) then +if ( ${fileName} == 'DIC_LOAD' || ${fileName} == 'DIC_VARS' ) then echo '#include "DIC_OPTIONS.h"' >> ${fileName}_temp endif if ( ${fileName} == 'GAD' ) then