Changeset 9208 for branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trcini.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/trcini.F90 (modified) (6 diffs)

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

@@ -28 +28 @@
    USE trcini_idtra    ! idealize tracer initialisation
    USE trcini_medusa   ! MEDUSA   initialisation
+   USE par_medusa      ! MEDUSA   parameters (needed for elemental cycles)
    USE trcdta          ! initialisation from files
    USE daymod          ! calendar manager
@@ -35 +36 @@
    USE sbc_oce
    USE trcice          ! tracers in sea ice
- 
+   USE sms_medusa      ! MEDUSA   initialisation
    IMPLICIT NONE
    PRIVATE
@@ -62 +63 @@
       !!                or read data or analytical formulation
       !!---------------------------------------------------------------------
-      INTEGER ::   jk, jn, jl    ! dummy loop indices
+      INTEGER ::   ji, jj, jk, jn, jl    ! dummy loop indices
+# if defined key_medusa && defined key_roam
+      !! AXY (23/11/2017)
+      REAL(wp)                         :: zsum3d, zsum2d
+      REAL(wp)                         :: zq1, zq2, loc_vol, loc_area
+      REAL(wp), DIMENSION(6)           :: loc_cycletot3, loc_cycletot2
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztot3d
+      REAL(wp), DIMENSION(jpi,jpj)     :: ztot2d, carea
+# endif
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -98 +107 @@
       !                                                              ! total volume of the ocean 
       areatot = glob_sum( cvol(:,:,:) )
+      carea(:,:) = e1e2t(:,:) * tmask(:,:,1) 
 
       IF( lk_pisces  )       CALL trc_ini_pisces       ! PISCES  bio-model
@@ -192 +202 @@
       ENDIF
 
+# if defined key_medusa && defined key_roam
+      ! AXY (17/11/2017): calculate initial totals of elemental cycles
+      !
+      ! This is done in a very hard-wired way here; in future, this could be
+      ! replaced with loops and using a 2D array; one dimension would cover
+      ! the tracers, the other would be for the elements; each tracer would
+      ! have a factor for each element to say how much of that element was
+      ! in that tracer; for example, PHN would be 1.0 for N, xrfn for Fe and
+      ! xthetapn for C, with the other elements 0.0; the array entry for PHN
+      ! would then be (1. 0. xrfn xthetapn 0. 0.) for (N, Si, Fe, C, A, O2);
+      ! saving this for the next iteration
+      !
+      cycletot(:) = 0._wp
+      ! report elemental totals at initialisation as we go along
+      IF ( lwp ) WRITE(numout,*)
+      IF ( lwp ) WRITE(numout,*)    ' Elemental cycle totals: '
+      IF ( lwp ) CALL flush(numout)
+      ! nitrogen
+      ztot3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
+                      trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin)
+      ztot2d(:,:)   = zn_sed_n(:,:)
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(1)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, cycletot(1)
+      IF ( lwp ) CALL flush(numout)
+      ! silicon
+      ztot3d(:,:,:) = trn(:,:,:,jppds) + trn(:,:,:,jpsil)
+      ztot2d(:,:)   = zn_sed_si(:,:)
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(2)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'silicon', zsum3d, zsum2d, cycletot(2)
+      IF ( lwp ) CALL flush(numout)
+      ! iron
+      ztot3d(:,:,:) = ((trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
+                      trn(:,:,:,jpzme) + trn(:,:,:,jpdet)) * xrfn) + trn(:,:,:,jpfer)
+      ztot2d(:,:)   = zn_sed_fe(:,:)
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(3)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'iron', zsum3d, zsum2d, cycletot(3)
+      IF ( lwp ) CALL flush(numout)
+      ! carbon (uses fixed C:N ratios on plankton tracers)
+      ztot3d(:,:,:) = (trn(:,:,:,jpphn) * xthetapn)  + (trn(:,:,:,jpphd) * xthetapd)  +  &
+                      (trn(:,:,:,jpzmi) * xthetazmi) + (trn(:,:,:,jpzme) * xthetazme) +  &
+                      trn(:,:,:,jpdtc) + trn(:,:,:,jpdic)
+      ztot2d(:,:)   = zn_sed_c(:,:) + zn_sed_ca(:,:)
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(4)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'carbon', zsum3d, zsum2d, cycletot(4)
+      IF ( lwp ) CALL flush(numout)
+      ! alkalinity (note benthic correction)
+      ztot3d(:,:,:) = trn(:,:,:,jpalk)
+      ztot2d(:,:)   = zn_sed_ca(:,:) * 2._wp
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(5)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'alkalinity', zsum3d, zsum2d, cycletot(5)
+      IF ( lwp ) CALL flush(numout)
+      ! oxygen (note no benthic)
+      ztot3d(:,:,:) = trn(:,:,:,jpoxy)
+      ztot2d(:,:)   = 0._wp
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(6)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'oxygen', zsum3d, zsum2d, cycletot(6)
+      IF ( lwp ) CALL flush(numout)
+      ! Check
+      zsum3d        = glob_sum( cvol(:,:,:) )
+      zsum2d        = glob_sum( carea(:,:) )      
+      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)
+      !
+# endif
+
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'trc_init : passive tracer set up completed'
@@ -202 +292 @@
 
 9000  FORMAT(' tracer nb : ',i2,'      name :',a10,'      initial content :',e18.10)
+9010  FORMAT(' element:',a10,                     &
+             ' 3d sum:',e18.10,' 2d sum:',e18.10, &
+             ' total:',e18.10)
       !
       IF( nn_timing == 1 )   CALL timing_stop('trc_init')