Changeset 9208 for branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trcrst.F90

Timestamp:

2018-01-10T17:06:13+01:00 (6 years ago)

Author:

jpalmier

Message:

JPALM -- GMED 364 -- clean/correct branch from GO6 Head

File:

• branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trcrst.F90 (modified) (4 diffs)

branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trcrst.F90

@@ -30 +30 @@
    USE daymod
    !! AXY (05/11/13): need these for MEDUSA to input/output benthic reservoirs
+   USE par_medusa
    USE sms_medusa
    USE trcsms_medusa
@@ -53 +54 @@
    PUBLIC   trc_rst_cal
    PUBLIC   trc_rst_stat
+#if defined key_medusa && defined key_roam
+   PUBLIC   trc_rst_conserve
+#endif
 
    !! * Substitutions
@@ -539 +543 @@
       IF( kt == nitrst ) THEN
           CALL trc_rst_stat            ! statistics
+#if defined key_medusa && defined key_roam
+          CALL trc_rst_conserve        ! conservation check
+#endif
           CALL iom_close( numrtw )     ! close the restart file (only at last time step)
 #if ! defined key_trdmxl_trc
@@ -705 +712 @@
 
 
+# if defined key_medusa && defined key_roam
+   SUBROUTINE trc_rst_conserve
+      !!----------------------------------------------------------------------
+      !!                    ***  trc_rst_conserve  ***
+      !!
+      !! ** purpose  :   Compute tracers conservation statistics
+      !!
+      !! AXY (17/11/2017)
+      !! This routine calculates the "now" inventories of the elemental 
+      !! cycles of MEDUSA and compares them to those calculate when the
+      !! model was initialised / restarted; the cycles calculated are:
+      !!    nitrogen, silicon, iron, carbon, alkalinity and oxygen
+      !!----------------------------------------------------------------------
+      INTEGER  :: ji, jj, jk, jn
+      REAL(wp) :: zsum3d, zsum2d, zinvt, zdelta, zratio, loc_vol, loc_are
+      REAL(wp) :: zq1, zq2, loc_vol, loc_area
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d, zvol
+      REAL(wp), DIMENSION(jpi,jpj)     :: z2d, zarea
+      REAL(wp), DIMENSION(6)           :: loc_cycletot3, loc_cycletot2
+      !!----------------------------------------------------------------------
+      !
+      IF( lwp ) THEN
+         WRITE(numout,*) 
+         WRITE(numout,*) '           ----TRACER CONSERVATION----             '
+         WRITE(numout,*) 
+      ENDIF
+      !
+      ! ocean volume
+      DO jk = 1, jpk
+         zvol(:,:,jk) = e1e2t(:,:) * fse3t_a(:,:,jk) * tmask(:,:,jk)
+      END DO
+      !
+      ! ocean area (for sediments)
+      zarea(:,:)      = e1e2t(:,:) * tmask(:,:,1)
+      !
+      !----------------------------------------------------------------------
+      ! nitrogen
+      z3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
+                   trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin)
+      z2d(:,:)   = zn_sed_n(:,:)
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(1)
+      zratio     = 1.0e2 * zdelta / cycletot(1)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, zinvt,   &
+         cycletot(1), zdelta, zratio
+      IF ( lwp ) WRITE(numout,*) 
+      !
+      !----------------------------------------------------------------------
+      ! silicon
+      z3d(:,:,:) = trn(:,:,:,jppds) + trn(:,:,:,jpsil)
+      z2d(:,:)   = zn_sed_si(:,:)
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(2)
+      zratio     = 1.0e2 * zdelta / cycletot(2)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'silicon', zsum3d, zsum2d, zinvt,    &
+         cycletot(2), zdelta, zratio
+      IF ( lwp ) WRITE(numout,*) 
+      !
+      !----------------------------------------------------------------------
+      ! iron
+      z3d(:,:,:) = ((trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
+            trn(:,:,:,jpzme) + trn(:,:,:,jpdet)) * xrfn) + trn(:,:,:,jpfer)
+      z2d(:,:)   = zn_sed_fe(:,:)
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(3)
+      zratio     = 1.0e2 * zdelta / cycletot(3)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'iron', zsum3d, zsum2d, zinvt,       &
+         cycletot(3), zdelta, zratio
+      IF ( lwp ) WRITE(numout,*) 
+      !
+      !----------------------------------------------------------------------
+      ! carbon
+      z3d(:,:,:) = (trn(:,:,:,jpphn) * xthetapn)  + (trn(:,:,:,jpphd) * xthetapd)  + &
+                   (trn(:,:,:,jpzmi) * xthetazmi) + (trn(:,:,:,jpzme) * xthetazme) + &
+                   trn(:,:,:,jpdtc) + trn(:,:,:,jpdic)
+      z2d(:,:)   = zn_sed_c(:,:) + zn_sed_ca(:,:)
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(4)
+      zratio     = 1.0e2 * zdelta / cycletot(4)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'carbon', zsum3d, zsum2d, zinvt,     &
+         cycletot(4), zdelta, zratio
+      IF ( lwp ) WRITE(numout,*) 
+      !
+      !----------------------------------------------------------------------
+      ! alkalinity
+      z3d(:,:,:) = trn(:,:,:,jpalk)
+      z2d(:,:)   = zn_sed_ca(:,:) * 2.0
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(5)
+      zratio     = 1.0e2 * zdelta / cycletot(5)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'alkalinity', zsum3d, zsum2d, zinvt, &
+         cycletot(5), zdelta, zratio
+      IF ( lwp ) WRITE(numout,*) 
+      !
+      !----------------------------------------------------------------------
+      ! oxygen
+      z3d(:,:,:) = trn(:,:,:,jpoxy)
+      z2d(:,:)   = 0.0
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(6)
+      zratio     = 1.0e2 * zdelta / cycletot(6)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'oxygen', zsum3d, zsum2d, zinvt,     &
+         cycletot(6), zdelta, zratio
+      !
+      !----------------------------------------------------------------------
+      ! Check 
+      zsum3d        = glob_sum( zvol(:,:,:) )
+      zsum2d        = glob_sum( zarea(:,:) )
+      IF ( lwp ) WRITE(numout,*)
+      IF ( lwp ) WRITE(numout,*) ' check : cvol    : ', zsum3d
+      IF ( lwp ) WRITE(numout,*) ' check : carea   : ', zsum2d
+      IF ( lwp ) WRITE(numout,*)
+      IF ( lwp ) CALL flush(numout)
+      !
+9010  FORMAT(' element:',a10,                     &
+             ' 3d sum:',e18.10,' 2d sum:',e18.10, &
+             ' total:',e18.10,' initial:',e18.10, &
+             ' delta:',e18.10,' %:',e18.10)
+      !
+   END SUBROUTINE trc_rst_conserve 
+# endif
+
+
 #else
    !!----------------------------------------------------------------------