changeset 89:85cec98f272f

Several bug fixes and improvements: + correct sinking speeds (per timestep) wpoco etc. for off-line (sediment only); + SED_OFFLINE_DEBUG saves extra bottom water output; + start off-line (sediment only) not before nday==2 of MICOM (otherwise last output is stored after spin-up); + rewind date more consistenly (and only for full MICOM/HAMOCC); + use INTENT(IN) for all *_ variables (and so use original vars for full); + sedfluxo may be updated also in off-line case; not sure if it makes a difference.
author Marco van Hulten <marco@hulten.org>
date Fri, 26 Jul 2019 16:03:00 +0200
parents 4d0b1c84ebea
children 918a4a892a21
files beleg_bgc.F90 mo_biomod.F90 mo_sedmnt_offline.F90 mod_dia.F powach.F90 sediment_step.F90
diffstat 6 files changed, 114 insertions(+), 31 deletions(-) [+]
line wrap: on
line diff
--- a/beleg_bgc.F90	Tue Jul 02 11:02:37 2019 +0200
+++ b/beleg_bgc.F90	Fri Jul 26 16:03:00 2019 +0200
@@ -191,10 +191,20 @@
       wpoc  =  5.*dtb       !m/d  iris : 5.
       wcal  = 30.*dtb       !m/d 
       wopal = 30.*dtb       !m/d  iris : 60
+# ifdef SED_OFFLINE
+      wpoco =  5.*dto       !m/d  iris : 5.
+      wcalo = 30.*dto       !m/d
+      wopalo= 30.*dto       !m/d  iris : 60
+# endif
 #ifdef WLIN
       wmin  =  1.*dtb       !m/d   minimum sinking speed
       wmax  = 60.*dtb       !m/d   maximum sinking speed
       wlin  = 60./2400.*dtb !m/d/m constant describing incr. with depth, r/a=1.0
+# ifdef SED_OFFLINE
+      wmino =  1.*dto       !m/d   minimum sinking speed
+      wmaxo = 60.*dto       !m/d   maximum sinking speed
+      wlino = 60./2400.*dto !m/d/m constant describing incr. with depth, r/a=1.0
+# endif
 #endif
 
       
@@ -410,6 +420,9 @@
      &         * (claydens - 1025.) / 1.567 * 1000.    &  !excess density / dyn. visc.
      &         * dustd2 * 1.e-4)*dtb
       wdust = dustsink
+#ifdef SED_OFFLINE
+      wdusto= dustsink * dto/dtb
+#endif
 
       IF (mnproc.eq.1) THEN
       WRITE(io_stdo_bgc,*)                                             &
--- a/mo_biomod.F90	Tue Jul 02 11:02:37 2019 +0200
+++ b/mo_biomod.F90	Fri Jul 26 16:03:00 2019 +0200
@@ -49,6 +49,9 @@
       REAL :: bluefix,tf2,tf1,tf0,tff  
       REAL :: bkphy,bkzoo,bkopal,bifr13,bifr14
       REAL :: wpoc,wcal,wopal,drempoc,dremdoc,dremn2o
+# if defined(SED_OFFLINE)
+      REAL :: wpoco,wcalo,wopalo,wdusto
+# endif
       REAL :: dphymor,dzoomor,dremopal
       REAL :: dremsul
       REAL :: psedi,csedi,ssedi
@@ -63,6 +66,9 @@
       REAL :: dustd1,dustd2,dustd3,dustsink,calmax
 #elif defined(WLIN)
       REAL :: wmin,wmax,wlin
+# if defined(SED_OFFLINE)
+      REAL :: wmino,wmaxo,wlino
+# endif
 #endif
 #ifdef __c_isotopes
       REAL :: factor_13c, factor_14c, atm_c14_cal, atm_dc14_cal
--- a/mo_sedmnt_offline.F90	Tue Jul 02 11:02:37 2019 +0200
+++ b/mo_sedmnt_offline.F90	Fri Jul 26 16:03:00 2019 +0200
@@ -75,8 +75,8 @@
 use mo_carbch,      only: ocetra, keqb, co3
 use mo_param1_bgc,  only: nocetra, idet, icalc, iopal, ifdust        &
    &                    , ks, nsedtra, npowtra
-use mo_biomod,      only: kbo, bolay, wpoc, wmin, wmax, wdust        &
-   &                    , wlin, wopal, wcal
+use mo_biomod,      only: kbo, bolay, wpoco, wmino, wmaxo, wdusto    &
+   &                    , wlino, wopalo, wcalo
 use mod_xc
 use mo_common_bgc
 use mod_dia,        only: ddm, depthslev, depthslev_bnds             &
@@ -291,11 +291,11 @@
       co3_kbo_avg = co3_kbo_avg / nstep_in_month
       do j = 1, jdm     ! Here we assume the WLIN case.
          do i = 1, idm
-            wpoc = min(wmin+wlin*bolay_avg(i,j),wmax)
-            prorca_avg(i,j)=ocetra_kbo_avg(i,j,idet  )*wpoc
-            prcaca_avg(i,j)=ocetra_kbo_avg(i,j,icalc )*wcal
-            silpro_avg(i,j)=ocetra_kbo_avg(i,j,iopal )*wopal
-            produs_avg(i,j)=ocetra_kbo_avg(i,j,ifdust)*wdust
+            wpoco = min(wmino+wlino*bolay_avg(i,j),wmaxo)
+            prorca_avg(i,j)=ocetra_kbo_avg(i,j,idet  )*wpoco
+            prcaca_avg(i,j)=ocetra_kbo_avg(i,j,icalc )*wcalo
+            silpro_avg(i,j)=ocetra_kbo_avg(i,j,iopal )*wopalo
+            produs_avg(i,j)=ocetra_kbo_avg(i,j,ifdust)*wdusto
          enddo
       enddo
       nstep_in_month = 0
@@ -391,6 +391,9 @@
    integer  :: nday_save, nmonth_save   ! calendar variables outside sedmnt_offline()
    integer  :: nyear_save, nday_of_year_save, nd_in_m_2, nday_in_year_save
    integer  :: nacc_bgc_1, nacc_bgc_2
+#if defined(SED_OFFLINE_DEBUG)
+   character(len= 2) seqstring_burst, seqstring_year
+#endif
 
    ! Read the bottom water climatology
    !
@@ -399,8 +402,11 @@
       lread_clim = .false.       !  but only one time.
    endif
 
-   if ( lsed_spinup .and. lcompleted_clim ) then
-      if (.not. (nmonth==1 .and. nday==1) .and. mnproc==1) then
+   ! Apparently we need to wait until nday = 2 before starting the off-line
+   ! spin-up, because otherwise data won't be saved for the last period.
+   !
+   if ( lsed_spinup .and. lcompleted_clim .and. nday==2 ) then
+      if (.not. (nmonth==1 .and. nday<=2) .and. mnproc==1) then
          write(io_stdo_bgc,*) 'sedmnt_offline(): WARNING: Not at start of year!'
       endif
 
@@ -435,6 +441,39 @@
                &         'sedmnt_offline(): nyear_global = ', nyear_global
 
          do nmonth = 1, 12
+#if defined(SED_OFFLINE_DEBUG)
+            ! We misuse the *_avg array variables to write current bottom water tracers.
+            ! An array of 12 timeslices is reserved, but only the first one stored (FIXME).
+            !
+            do iocetra = 1, nocetra
+               do j = 1, jdm
+                  do i = 1, idm
+                     ocetra_kbo_avg(i,j,iocetra) = ocetra(i,j,kbo(i,j),iocetra)
+                  enddo
+               enddo
+            enddo
+            do j = 1, jdm
+               do i = 1, idm
+                  bgc_t_kbo_avg(i,j) = bgc_t(i,j,kbo(i,j))
+                  bgc_s_kbo_avg(i,j) = bgc_s(i,j,kbo(i,j))
+                  bgc_rho_kbo_avg(i,j) = bgc_rho(i,j,kbo(i,j))
+                  co3_kbo_avg(i,j) = co3(i,j,kbo(i,j))
+                  wpoco = min(wmino+wlino*bolay_avg(i,j),wmaxo)
+                  prorca_avg(i,j)=ocetra_kbo_avg(i,j,idet  )*wpoco
+                  prcaca_avg(i,j)=ocetra_kbo_avg(i,j,icalc )*wcalo
+                  silpro_avg(i,j)=ocetra_kbo_avg(i,j,iopal )*wopalo
+                  produs_avg(i,j)=ocetra_kbo_avg(i,j,ifdust)*wdusto
+               enddo
+            enddo
+            bolay_avg = bolay
+            keqb_avg  = keqb
+
+            write (seqstring_burst,'(I0.2)') nburst
+            write (seqstring_year,'(I0.2)') nyear
+            call aufw_bgc_onlysed(idm, jdm, kdm, nyear, nmonth, nday, nstep,        &
+               & "bottom_seawater.burst"//seqstring_burst//".year"//seqstring_year//".DEBUG.nc")
+#endif
+
             ! set up sediment layers (much higher diffusion rate; correct timestep)
             dtoff = 3600*24*nd_in_m(nmonth)
             call bodensed(kpie,kpje,kpke,pddpo)
@@ -503,8 +542,24 @@
       powtra2(:,:,ks+1:2*ks,:) = powtra(:,:,:,:)
       burial2(:,:,1,:)         = burial(:,:,:)
       burial2(:,:,2,:)         = burial(:,:,:)
+
+#if defined(SED_OFFLINE_DEBUG)
+      ! We misused the *_avg array variables to write current bottom water tracers.
+      !
+      ocetra_kbo_avg = 0.0
+      bolay_avg = 4000.0
+      keqb_avg  = 0.0
+      prorca_avg = 0.0
+      prcaca_avg = 0.0
+      silpro_avg = 0.0
+      produs_avg = 0.0
+      bgc_t_kbo_avg = 0.0
+      bgc_s_kbo_avg = 0.0
+      bgc_rho_kbo_avg = 0.0
+      co3_kbo_avg = 0.0
+#endif
+      lcompleted_clim = .false.  ! don't reuse the climatology
    endif
-   lcompleted_clim = .false.  ! don't reuse the climatology
 end subroutine sedmnt_offline
 
 subroutine alloc_mem_sedmnt_offline(kpie, kpje)
--- a/mod_dia.F	Tue Jul 02 11:02:37 2019 +0200
+++ b/mod_dia.F	Fri Jul 26 16:03:00 2019 +0200
@@ -228,6 +228,7 @@
       integer runid_len,nstep
       real diagfq,filefq
       logical filemon,fileann,filedec,filecen,filemil
+      integer mymonth
 c
       real epsil
       parameter (epsil=1.e-11)
@@ -265,10 +266,18 @@
      .        runid(1:runid_len),'_',ctag,'_',ny,'.',nm,'.',nd,'.nc'
           endif
         elseif (expcnf.eq.'cesm') then
-          if ( filemon .or. fileann ) then
-            call pstdat(ny,nm,nd,1)
-          else
-            call pstdat(ny,nm,nd,nint(diagfq))
+          if (.not. lspinning_up_sed ) then
+            if (diagfq==30) then ! always rewind the actual diagfq as the
+                                 !  fname is created after the first diagfq
+              if (nmonth /= 1) then
+                 mymonth = nmonth-1
+              else
+                 mymonth = 12
+              endif
+              call pstdat(ny,nm,nd,nd_in_m(mymonth))
+            else
+              call pstdat(ny,nm,nd,nint(diagfq))
+            endif
           endif
           if (filemon) then
             if ( lspinning_up_sed ) then
--- a/powach.F90	Tue Jul 02 11:02:37 2019 +0200
+++ b/powach.F90	Fri Jul 26 16:03:00 2019 +0200
@@ -47,10 +47,10 @@
 !  REAL     bolay    - thickness of the bottom gridbox [m].
 !  REAL     ocetra   - bottom layer ocean tracer ocetra(:,:,kbo,:).
 !  REAL     keqb     - chemical equilibrium constants.
-!  REAL     prorca   - sedimentation of carbon [kmol/m^2/s].
-!  REAL     prcaca   - sedimentation of calcium carbonate [kmol/m^2/s].
-!  REAL     silpro   - sedimentation of biogenic silica [kmol/m^2/s].
-!  REAL     produs   - sedimentation of lithogenic dust [kmol/m^2/s].
+!  REAL     prorca   - sedimentation of carbon [kmol/m^2/timestep].
+!  REAL     prcaca   - sedimentation of calcium carbonate [kmol/m^2/timestep].
+!  REAL     silpro   - sedimentation of biogenic silica [kmol/m^2/timestep].
+!  REAL     produs   - sedimentation of lithogenic dust [kmol/m^2/timestep].
 !  REAL     co3      - dissolved carbonate in the bottom gridbox [mol/kg].
 !
 ! The actual argument names sometimes end with an underscore, signifying
@@ -86,10 +86,10 @@
 real, intent(in)     :: bolay_ (kpie,kpje)
 real, intent(in)     :: ocetra_(kpie,kpje,nocetra)
 real, intent(in)     :: keqb_  (11,kpie,kpje)
-real, intent(inout)  :: prorca_(kpie,kpje)
-real, intent(inout)  :: prcaca_(kpie,kpje)
-real, intent(inout)  :: silpro_(kpie,kpje)
-real, intent(inout)  :: produs_(kpie,kpje)
+real, intent(in)     :: prorca_(kpie,kpje)
+real, intent(in)     :: prcaca_(kpie,kpje)
+real, intent(in)     :: silpro_(kpie,kpje)
+real, intent(in)     :: produs_(kpie,kpje)
 real, intent(in)     :: co3_   (kpie,kpje)
 real, intent(in)     :: omask  (kpie,kpje)
 
@@ -201,10 +201,10 @@
 
 do i = 1, kpie
    if(omask(i,j) > 0.5) then
-      if ( .not. lspinning_up_sed ) then
          sedfluxo(i,j,ipowasi) = sedfluxo(i,j,ipowasi) +                   &
             &    (silsat-sediso(i,0) - ocetra_(i,j,isilica))*bolay_(i,j)
 
+      if ( .not. lspinning_up_sed ) then
          ! NOTE: If ocetra_(:,:,:) instead of ocetra(:,:,kbo,:) were to be updated,
          !       the latter would need to be updated by the former later on (MvH)!
          !
@@ -562,16 +562,16 @@
 !$OMP PARALLEL DO
 do j = 1, kpje
    do i = 1, kpie
-      silpro_(i,j) = 0.
-      prorca_(i,j) = 0.
+      silpro(i,j) = 0.
+      prorca(i,j) = 0.
 #ifdef __c_isotopes
       pror13(i,j) = 0.
       pror14(i,j) = 0.
       prca13(i,j) = 0.
       prca14(i,j) = 0.
 #endif
-      prcaca_(i,j) = 0.
-      produs_(i,j) = 0.
+      prcaca(i,j) = 0.
+      produs(i,j) = 0.
    enddo
 enddo
 !$OMP END PARALLEL DO
--- a/sediment_step.F90	Tue Jul 02 11:02:37 2019 +0200
+++ b/sediment_step.F90	Fri Jul 26 16:03:00 2019 +0200
@@ -71,10 +71,10 @@
 real, intent(in)     :: ocetra_(kpie,kpje,nocetra)
 real, intent(in)     :: bolay_ (kpie,kpje)
 real, intent(in)     :: keqb_  (11,kpie,kpje)
-real, intent(inout)  :: prorca_(kpie,kpje)
-real, intent(inout)  :: prcaca_(kpie,kpje)
-real, intent(inout)  :: silpro_(kpie,kpje)
-real, intent(inout)  :: produs_(kpie,kpje)
+real, intent(in)     :: prorca_(kpie,kpje)
+real, intent(in)     :: prcaca_(kpie,kpje)
+real, intent(in)     :: silpro_(kpie,kpje)
+real, intent(in)     :: produs_(kpie,kpje)
 real, intent(in)     :: co3_   (kpie,kpje)
 integer, intent(in), optional :: jpowaic_(nbgcmax)
 integer, intent(in), optional :: jpowaal_(nbgcmax)