Changeset 8926 for branches/NERC


Ignore:
Timestamp:
2017-12-06T17:10:46+01:00 (4 years ago)
Author:
jpalmier
Message:

JPALM — add conservation checks inside MEDUSA

Location:
branches/NERC/dev_r5518_GO6_conserv_Check/NEMOGCM/NEMO
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r5518_GO6_conserv_Check/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r6823 r8926  
    851851      CALL wrk_alloc( jpi,jpj, rgt33 ) 
    852852      ! 
     853#if defined key_coare_bulk_patch 
     854     rgt33 = MIN(zw10, 22.) 
     855     !! cd_neutral_10m = 1.e-3 * ( 2.7/rgt33  + 0.142 + 0.06*rgt33  + 
     856     !!                  0.0025*rgt33 **2 - 1.25e-9*rgt33 **6) 
     857     !! A new version May 31st! 
     858      cd_neutral_10m = 1.e-3 * ( 2.4/rgt33  + 0.15 + 0.095*rgt33  + 0.0008*rgt33 **2 - 1.0e-9*rgt33 **6) 
     859#else 
    853860      !! When wind speed > 33 m/s => Cyclone conditions => special treatment 
    854861      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1   
     
    856863         &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. 
    857864         &      + rgt33         *      2.34   )                                                    ! zw10 >= 33. 
     865#endif 
    858866      ! 
    859867      CALL wrk_dealloc( jpi,jpj, rgt33) 
  • branches/NERC/dev_r5518_GO6_conserv_Check/NEMOGCM/NEMO/TOP_SRC/trc.F90

    r8521 r8926  
    3333   !! -------------------------------------------------- 
    3434   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)               ::  trai           !: initial total tracer 
     35# if defined key_medusa && key_roam  
     36   !! AXY (17/11/2017): elemental cycle initial totals 
     37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)               ::  cycletot       !: initial elemental cycle total 
     38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)               ::  cycletot2      !: initial elemental cycle total excl. halo in mpp_sum 
     39# endif 
    3540   REAL(wp), PUBLIC                                                ::  areatot        !: total volume  
    3641   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:  )         ::  cvol           !: volume correction -degrad option-  
     
    266271         &      cvol(jpi,jpj,jpk)     , rdttrc(jpk)           , trai(jptra)           ,       & 
    267272         &      ctrcnm(jptra)         , ctrcln(jptra)         , ctrcun(jptra)         ,       &  
     273# if defined key_medusa && defined key_roam 
     274         &      cycletot(6), cycletot2(6)                                             ,       & 
     275# endif 
    268276         &      ln_trc_ini(jptra)     , ln_trc_wri(jptra)     , qsr_mean(jpi,jpj)     ,  STAT = trc_alloc  )   
    269277 
  • branches/NERC/dev_r5518_GO6_conserv_Check/NEMOGCM/NEMO/TOP_SRC/trcini.F90

    r8442 r8926  
    2828   USE trcini_idtra    ! idealize tracer initialisation 
    2929   USE trcini_medusa   ! MEDUSA   initialisation 
     30   USE par_medusa      ! MEDUSA   parameters (needed for elemental cycles) 
    3031   USE trcdta          ! initialisation from files 
    3132   USE daymod          ! calendar manager 
     
    3536   USE sbc_oce 
    3637   USE trcice          ! tracers in sea ice 
    37   
     38   USE sms_medusa      ! MEDUSA   initialisation 
    3839   IMPLICIT NONE 
    3940   PRIVATE 
     
    6263      !!                or read data or analytical formulation 
    6364      !!--------------------------------------------------------------------- 
    64       INTEGER ::   jk, jn, jl    ! dummy loop indices 
     65      INTEGER ::   ji, jj, jk, jn, jl    ! dummy loop indices 
     66# if defined key_medusa && defined key_roam 
     67      !! AXY (23/11/2017) 
     68      REAL(wp)                         :: zsum3d, zsum2d 
     69      REAL(wp)                         :: zq1, zq2, loc_vol, loc_area 
     70      REAL(wp), DIMENSION(6)           :: loc_cycletot3, loc_cycletot2 
     71      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztot3d 
     72      REAL(wp), DIMENSION(jpi,jpj)     :: ztot2d, carea 
     73# endif 
    6574      CHARACTER (len=25) :: charout 
    6675      !!--------------------------------------------------------------------- 
     
    98107      !                                                              ! total volume of the ocean  
    99108      areatot = glob_sum( cvol(:,:,:) ) 
     109      carea(:,:) = e1e2t(:,:) * tmask(:,:,1)  
    100110 
    101111      IF( lk_pisces  )       CALL trc_ini_pisces       ! PISCES  bio-model 
     
    192202      ENDIF 
    193203 
     204# if defined key_medusa && defined key_roam 
     205      ! AXY (17/11/2017): calculate initial totals of elemental cycles 
     206      ! 
     207      ! This is done in a very hard-wired way here; in future, this could be 
     208      ! replaced with loops and using a 2D array; one dimension would cover 
     209      ! the tracers, the other would be for the elements; each tracer would 
     210      ! have a factor for each element to say how much of that element was 
     211      ! in that tracer; for example, PHN would be 1.0 for N, xrfn for Fe and 
     212      ! xthetapn for C, with the other elements 0.0; the array entry for PHN 
     213      ! would then be (1. 0. xrfn xthetapn 0. 0.) for (N, Si, Fe, C, A, O2); 
     214      ! saving this for the next iteration 
     215      ! 
     216      cycletot(:) = 0._wp 
     217      ! report elemental totals at initialisation as we go along 
     218      IF ( lwp ) WRITE(numout,*) 
     219      IF ( lwp ) WRITE(numout,*)    ' Elemental cycle totals: ' 
     220      IF ( lwp ) CALL flush(numout) 
     221      ! nitrogen 
     222      ztot3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 
     223                      trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin) 
     224      ztot2d(:,:)   = zn_sed_n(:,:) 
     225      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 
     226      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) ) 
     227      cycletot(1)   = zsum3d + zsum2d 
     228      IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, cycletot(1) 
     229      IF ( lwp ) CALL flush(numout) 
     230      ! silicon 
     231      ztot3d(:,:,:) = trn(:,:,:,jppds) + trn(:,:,:,jpsil) 
     232      ztot2d(:,:)   = zn_sed_si(:,:) 
     233      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 
     234      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) ) 
     235      cycletot(2)   = zsum3d + zsum2d 
     236      IF ( lwp ) WRITE(numout,9010) 'silicon', zsum3d, zsum2d, cycletot(2) 
     237      IF ( lwp ) CALL flush(numout) 
     238      ! iron 
     239      ztot3d(:,:,:) = ((trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 
     240                      trn(:,:,:,jpzme) + trn(:,:,:,jpdet)) * xrfn) + trn(:,:,:,jpfer) 
     241      ztot2d(:,:)   = zn_sed_fe(:,:) 
     242      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 
     243      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) ) 
     244      cycletot(3)   = zsum3d + zsum2d 
     245      IF ( lwp ) WRITE(numout,9010) 'iron', zsum3d, zsum2d, cycletot(3) 
     246      IF ( lwp ) CALL flush(numout) 
     247      ! carbon (uses fixed C:N ratios on plankton tracers) 
     248      ztot3d(:,:,:) = (trn(:,:,:,jpphn) * xthetapn)  + (trn(:,:,:,jpphd) * xthetapd)  +  & 
     249                      (trn(:,:,:,jpzmi) * xthetazmi) + (trn(:,:,:,jpzme) * xthetazme) +  & 
     250                      trn(:,:,:,jpdtc) + trn(:,:,:,jpdic) 
     251      ztot2d(:,:)   = zn_sed_c(:,:) + zn_sed_ca(:,:) 
     252      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 
     253      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) ) 
     254      cycletot(4)   = zsum3d + zsum2d 
     255      IF ( lwp ) WRITE(numout,9010) 'carbon', zsum3d, zsum2d, cycletot(4) 
     256      IF ( lwp ) CALL flush(numout) 
     257      ! alkalinity (note benthic correction) 
     258      ztot3d(:,:,:) = trn(:,:,:,jpalk) 
     259      ztot2d(:,:)   = zn_sed_ca(:,:) * 2._wp 
     260      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 
     261      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) ) 
     262      cycletot(5)   = zsum3d + zsum2d 
     263      IF ( lwp ) WRITE(numout,9010) 'alkalinity', zsum3d, zsum2d, cycletot(5) 
     264      IF ( lwp ) CALL flush(numout) 
     265      ! oxygen (note no benthic) 
     266      ztot3d(:,:,:) = trn(:,:,:,jpoxy) 
     267      ztot2d(:,:)   = 0._wp 
     268      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) ) 
     269      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) ) 
     270      cycletot(6)   = zsum3d + zsum2d 
     271      IF ( lwp ) WRITE(numout,9010) 'oxygen', zsum3d, zsum2d, cycletot(6) 
     272      IF ( lwp ) CALL flush(numout) 
     273      ! Check 
     274      zsum3d        = glob_sum( cvol(:,:,:) ) 
     275      zsum2d        = glob_sum( carea(:,:) )       
     276      IF ( lwp ) WRITE(numout,*) 
     277      IF ( lwp ) WRITE(numout,*) ' check : cvol    : ', zsum3d 
     278      IF ( lwp ) WRITE(numout,*) ' check : carea   : ', zsum2d 
     279      IF ( lwp ) WRITE(numout,*) 
     280      IF ( lwp ) CALL flush(numout) 
     281      ! 
     282      ! This second check dodges around the halo points in the grid 
     283      ! to check that glob_sum is doing what it's supposed to be 
     284      ! doing; note that the loop ordering here is the "correct" way  
     285      ! (according to Dr. Google) 
     286      ! 
     287      IF ( lwp ) WRITE(numout,*) 
     288      IF ( lwp ) WRITE(numout,*)    ' Elemental cycle totals (check): ' 
     289      ! ocean cells 
     290      loc_cycletot3(:) = 0._wp 
     291      loc_vol          = 0._wp 
     292      DO jk = 1, jpk 
     293         DO jj = nldj,nlej 
     294            DO ji = nldi,nlei 
     295               loc_vol = loc_vol + cvol(ji,jj,jk) 
     296               ! nitrogen 
     297               zq1 = trn(ji,jj,jk,jpphn) + trn(ji,jj,jk,jpphd) + trn(ji,jj,jk,jpzmi) + & 
     298                     trn(ji,jj,jk,jpzme) + trn(ji,jj,jk,jpdet) + trn(ji,jj,jk,jpdin) 
     299               loc_cycletot3(1) = loc_cycletot3(1) + ( zq1 * cvol(ji,jj,jk) ) 
     300               ! silicon 
     301               zq1 = trn(ji,jj,jk,jppds) + trn(ji,jj,jk,jpsil) 
     302               loc_cycletot3(2) = loc_cycletot3(2) + ( zq1 * cvol(ji,jj,jk) ) 
     303               ! iron 
     304               zq1 = ((trn(ji,jj,jk,jpphn) + trn(ji,jj,jk,jpphd) + trn(ji,jj,jk,jpzmi) + & 
     305                     trn(ji,jj,jk,jpzme) + trn(ji,jj,jk,jpdet)) * xrfn) + trn(ji,jj,jk,jpfer) 
     306               loc_cycletot3(3) = loc_cycletot3(3) + ( zq1 * cvol(ji,jj,jk) ) 
     307               ! carbon 
     308               zq1 = (trn(ji,jj,jk,jpphn) * xthetapn)  + (trn(ji,jj,jk,jpphd) * xthetapd)  +  & 
     309                     (trn(ji,jj,jk,jpzmi) * xthetazmi) + (trn(ji,jj,jk,jpzme) * xthetazme) +  & 
     310                     trn(ji,jj,jk,jpdtc) + trn(ji,jj,jk,jpdic) 
     311               loc_cycletot3(4) = loc_cycletot3(4) + ( zq1 * cvol(ji,jj,jk) ) 
     312               ! alkalinity 
     313               zq1 = trn(ji,jj,jk,jpalk) 
     314               loc_cycletot3(5) = loc_cycletot3(5) + ( zq1 * cvol(ji,jj,jk) ) 
     315               ! oxygen 
     316               zq1 = trn(ji,jj,jk,jpoxy) 
     317               loc_cycletot3(6) = loc_cycletot3(6) + ( zq1 * cvol(ji,jj,jk) ) 
     318            ENDDO 
     319         ENDDO 
     320      ENDDO 
     321      ! 
     322      ! sediment cells 
     323      loc_cycletot2(:) = 0._wp 
     324      loc_area         = 0._wp 
     325      jk = 1 
     326      DO jj = nldj,nlej 
     327         DO ji = nldi,nlei 
     328            loc_area = loc_area + (e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     329            ! nitrogen 
     330            zq1 = zn_sed_n(ji,jj) 
     331            loc_cycletot2(1) = loc_cycletot2(1) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     332            ! silicon 
     333            zq1 = zn_sed_si(ji,jj) 
     334            loc_cycletot2(2) = loc_cycletot2(2) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     335            ! iron 
     336            zq1 = zn_sed_fe(ji,jj) 
     337            loc_cycletot2(3) = loc_cycletot2(3) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     338            ! carbon 
     339            zq1 = zn_sed_c(ji,jj) + zn_sed_ca(ji,jj) 
     340            loc_cycletot2(4) = loc_cycletot2(4) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     341            ! alkalinity 
     342            zq1 = zn_sed_ca(ji,jj) * 2._wp 
     343            loc_cycletot2(5) = loc_cycletot2(5) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     344            ! skip oxygen 
     345         ENDDO 
     346      ENDDO 
     347      ! 
     348      ! report back 
     349      zq1 = loc_cycletot3(1) 
     350      zq2 = loc_cycletot2(1)  
     351      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     352      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     353      cycletot2(1) = zq1 + zq2 
     354      IF ( lwp ) WRITE(numout,9010) 'nitrogen',   zq1, zq2, cycletot2(1) 
     355      zq1 = loc_cycletot3(2)  
     356      zq2 = loc_cycletot2(2)  
     357      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     358      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     359      cycletot2(2) = zq1 + zq2 
     360      IF ( lwp ) WRITE(numout,9010) 'silicon',    zq1, zq2, cycletot2(2) 
     361      zq1 = loc_cycletot3(3)  
     362      zq2 = loc_cycletot2(3)  
     363      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     364      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     365      cycletot2(3) = zq1 + zq2 
     366      IF ( lwp ) WRITE(numout,9010) 'iron',       zq1, zq2, cycletot2(3) 
     367      zq1 = loc_cycletot3(4)  
     368      zq2 = loc_cycletot2(4)  
     369      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     370      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     371      cycletot2(4) = zq1 + zq2 
     372      IF ( lwp ) WRITE(numout,9010) 'carbon',     zq1, zq2, cycletot2(4) 
     373      zq1 = loc_cycletot3(5)  
     374      zq2 = loc_cycletot2(5)  
     375      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     376      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     377      cycletot2(5) = zq1 + zq2 
     378      IF ( lwp ) WRITE(numout,9010) 'alkalinity', zq1, zq2, cycletot2(5) 
     379      zq1 = loc_cycletot3(6)  
     380      zq2 = loc_cycletot2(6)  
     381      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     382      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     383      cycletot2(6) = zq1 + zq2 
     384      IF ( lwp ) WRITE(numout,9010) 'oxygen',     zq1, zq2, cycletot2(6) 
     385      zq1 = loc_vol  
     386      zq2 = loc_area        
     387      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     388      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     389      IF ( lwp ) WRITE(numout,*) 
     390      IF ( lwp ) WRITE(numout,*) ' check : cvol    : ', zq1 
     391      IF ( lwp ) WRITE(numout,*) ' check : carea   : ', zq2 
     392      IF ( lwp ) WRITE(numout,*) 
     393      IF ( lwp ) CALL flush(numout) 
     394      ! 
     395      ! report back 
     396      ! nitrogen 
     397      ztot3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 
     398                      trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin) 
     399      ztot2d(:,:)   = zn_sed_n(:,:) 
     400      zsum3d        = glob_sum( ztot3d(nldi:nlei,nldj:nlej,:) * cvol(nldi:nlei,nldj:nlej,:) ) 
     401      zsum2d        = glob_sum( ztot2d(nldi:nlei,nldj:nlej) * carea(nldi:nlei,nldj:nlej) ) 
     402      cycletot(1)   = zsum3d + zsum2d 
     403      IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, cycletot(1) 
     404      IF ( lwp ) CALL flush(numout) 
     405 
     406      ! 
     407 
     408# endif 
     409 
    194410      IF(lwp) WRITE(numout,*) 
    195411      IF(lwp) WRITE(numout,*) 'trc_init : passive tracer set up completed' 
     
    202418 
    2034199000  FORMAT(' tracer nb : ',i2,'      name :',a10,'      initial content :',e18.10) 
     4209010  FORMAT(' element:',a10,                     & 
     421             ' 3d sum:',e18.10,' 2d sum:',e18.10, & 
     422             ' total:',e18.10) 
    204423      ! 
    205424      IF( nn_timing == 1 )   CALL timing_stop('trc_init') 
  • branches/NERC/dev_r5518_GO6_conserv_Check/NEMOGCM/NEMO/TOP_SRC/trcrst.F90

    r8427 r8926  
    3030   USE daymod 
    3131   !! AXY (05/11/13): need these for MEDUSA to input/output benthic reservoirs 
     32   USE par_medusa 
    3233   USE sms_medusa 
    3334   USE trcsms_medusa 
     
    5253   PUBLIC   trc_rst_cal 
    5354   PUBLIC   trc_rst_stat 
     55#if defined key_medusa && defined key_roam 
     56   PUBLIC   trc_rst_conserve 
     57#endif 
    5458   PUBLIC   trc_rst_dia_stat 
    5559   PUBLIC   trc_rst_tra_stat 
     
    531535      IF( kt == nitrst ) THEN 
    532536          CALL trc_rst_stat            ! statistics 
     537#if defined key_medusa && defined key_roam 
     538          CALL trc_rst_conserve        ! conservation check 
     539#endif 
    533540          CALL iom_close( numrtw )     ! close the restart file (only at last time step) 
    534541#if ! defined key_trdmxl_trc 
     
    697704 
    698705 
     706#if defined key_medusa && defined key_roam 
     707   SUBROUTINE trc_rst_conserve 
     708      !!---------------------------------------------------------------------- 
     709      !!                    ***  trc_rst_conserve  *** 
     710      !! 
     711      !! ** purpose  :   Compute tracers conservation statistics 
     712      !! 
     713      !! AXY (17/11/2017) 
     714      !! This routine calculates the "now" inventories of the elemental  
     715      !! cycles of MEDUSA and compares them to those calculate when the 
     716      !! model was initialised / restarted; the cycles calculated are: 
     717      !!    nitrogen, silicon, iron, carbon, alkalinity and oxygen 
     718      !!---------------------------------------------------------------------- 
     719      INTEGER  :: ji, jj, jk, jn 
     720      REAL(wp) :: zsum3d, zsum2d, zinvt, zdelta, zratio, loc_vol, loc_are 
     721      REAL(wp) :: zq1, zq2, loc_vol, loc_area 
     722      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d, zvol 
     723      REAL(wp), DIMENSION(jpi,jpj)     :: z2d, zarea 
     724      REAL(wp), DIMENSION(6)           :: loc_cycletot3, loc_cycletot2 
     725      !!---------------------------------------------------------------------- 
     726      ! 
     727      IF( lwp ) THEN 
     728         WRITE(numout,*)  
     729         WRITE(numout,*) '           ----TRACER CONSERVATION----             ' 
     730         WRITE(numout,*)  
     731      ENDIF 
     732      ! 
     733      ! ocean volume 
     734      DO jk = 1, jpk 
     735         zvol(:,:,jk) = e1e2t(:,:) * fse3t_a(:,:,jk) * tmask(:,:,jk) 
     736      END DO 
     737      ! 
     738      ! ocean area (for sediments) 
     739      zarea(:,:)      = e1e2t(:,:) * tmask(:,:,1) 
     740      ! 
     741      !---------------------------------------------------------------------- 
     742      ! nitrogen 
     743      z3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 
     744                   trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin) 
     745      z2d(:,:)   = zn_sed_n(:,:) 
     746      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 
     747      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) ) 
     748      ! total tracer, and delta 
     749      zinvt      = zsum3d + zsum2d 
     750      zdelta     = zinvt - cycletot(1) 
     751      zratio     = 1.0e2 * zdelta / cycletot(1) 
     752      ! 
     753      IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, zinvt,   & 
     754         cycletot(1), zdelta, zratio 
     755      IF ( lwp ) WRITE(numout,*)  
     756      ! 
     757      !---------------------------------------------------------------------- 
     758      ! silicon 
     759      z3d(:,:,:) = trn(:,:,:,jppds) + trn(:,:,:,jpsil) 
     760      z2d(:,:)   = zn_sed_si(:,:) 
     761      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 
     762      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) ) 
     763      ! total tracer, and delta 
     764      zinvt      = zsum3d + zsum2d 
     765      zdelta     = zinvt - cycletot(2) 
     766      zratio     = 1.0e2 * zdelta / cycletot(2) 
     767      ! 
     768      IF ( lwp ) WRITE(numout,9010) 'silicon', zsum3d, zsum2d, zinvt,    & 
     769         cycletot(2), zdelta, zratio 
     770      IF ( lwp ) WRITE(numout,*)  
     771      ! 
     772      !---------------------------------------------------------------------- 
     773      ! iron 
     774      z3d(:,:,:) = ((trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 
     775            trn(:,:,:,jpzme) + trn(:,:,:,jpdet)) * xrfn) + trn(:,:,:,jpfer) 
     776      z2d(:,:)   = zn_sed_fe(:,:) 
     777      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 
     778      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) ) 
     779      ! total tracer, and delta 
     780      zinvt      = zsum3d + zsum2d 
     781      zdelta     = zinvt - cycletot(3) 
     782      zratio     = 1.0e2 * zdelta / cycletot(3) 
     783      ! 
     784      IF ( lwp ) WRITE(numout,9010) 'iron', zsum3d, zsum2d, zinvt,       & 
     785         cycletot(3), zdelta, zratio 
     786      IF ( lwp ) WRITE(numout,*)  
     787      ! 
     788      !---------------------------------------------------------------------- 
     789      ! carbon 
     790      z3d(:,:,:) = (trn(:,:,:,jpphn) * xthetapn)  + (trn(:,:,:,jpphd) * xthetapd)  + & 
     791                   (trn(:,:,:,jpzmi) * xthetazmi) + (trn(:,:,:,jpzme) * xthetazme) + & 
     792                   trn(:,:,:,jpdtc) + trn(:,:,:,jpdic) 
     793      z2d(:,:)   = zn_sed_c(:,:) + zn_sed_ca(:,:) 
     794      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 
     795      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) ) 
     796      ! total tracer, and delta 
     797      zinvt      = zsum3d + zsum2d 
     798      zdelta     = zinvt - cycletot(4) 
     799      zratio     = 1.0e2 * zdelta / cycletot(4) 
     800      ! 
     801      IF ( lwp ) WRITE(numout,9010) 'carbon', zsum3d, zsum2d, zinvt,     & 
     802         cycletot(4), zdelta, zratio 
     803      IF ( lwp ) WRITE(numout,*)  
     804      ! 
     805      !---------------------------------------------------------------------- 
     806      ! alkalinity 
     807      z3d(:,:,:) = trn(:,:,:,jpalk) 
     808      z2d(:,:)   = zn_sed_ca(:,:) * 2.0 
     809      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 
     810      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) ) 
     811      ! total tracer, and delta 
     812      zinvt      = zsum3d + zsum2d 
     813      zdelta     = zinvt - cycletot(5) 
     814      zratio     = 1.0e2 * zdelta / cycletot(5) 
     815      ! 
     816      IF ( lwp ) WRITE(numout,9010) 'alkalinity', zsum3d, zsum2d, zinvt, & 
     817         cycletot(5), zdelta, zratio 
     818      IF ( lwp ) WRITE(numout,*)  
     819      ! 
     820      !---------------------------------------------------------------------- 
     821      ! oxygen 
     822      z3d(:,:,:) = trn(:,:,:,jpoxy) 
     823      z2d(:,:)   = 0.0 
     824      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 
     825      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) ) 
     826      ! total tracer, and delta 
     827      zinvt      = zsum3d + zsum2d 
     828      zdelta     = zinvt - cycletot(6) 
     829      zratio     = 1.0e2 * zdelta / cycletot(6) 
     830      ! 
     831      IF ( lwp ) WRITE(numout,9010) 'oxygen', zsum3d, zsum2d, zinvt,     & 
     832         cycletot(6), zdelta, zratio 
     833      ! 
     834      !---------------------------------------------------------------------- 
     835      ! Check  
     836      zsum3d        = glob_sum( zvol(:,:,:) ) 
     837      zsum2d        = glob_sum( zarea(:,:) ) 
     838      IF ( lwp ) WRITE(numout,*) 
     839      IF ( lwp ) WRITE(numout,*) ' check : cvol    : ', zsum3d 
     840      IF ( lwp ) WRITE(numout,*) ' check : carea   : ', zsum2d 
     841      IF ( lwp ) WRITE(numout,*) 
     842      IF ( lwp ) CALL flush(numout) 
     843      ! 
     844      !---------------------------------------------------------------------- 
     845      ! Excluding Halo version: 
     846      ! 
     847      ! This second check dodges around the halo points in the grid 
     848      ! to check that glob_sum is doing what it's supposed to be 
     849      ! doing; note that the loop ordering here is the "correct" way  
     850      ! (according to Dr. Google) 
     851      ! 
     852      IF ( lwp ) WRITE(numout,*) 
     853      IF ( lwp ) WRITE(numout,*)    ' Elemental cycle totals (check): ' 
     854      ! ocean cells 
     855      loc_cycletot3(:) = 0._wp 
     856      loc_vol          = 0._wp 
     857      DO jk = 1, jpk 
     858         DO jj = nldj,nlej 
     859            DO ji = nldi,nlei 
     860               loc_vol = loc_vol + cvol(ji,jj,jk) 
     861               ! nitrogen 
     862               zq1 = trn(ji,jj,jk,jpphn) + trn(ji,jj,jk,jpphd) + trn(ji,jj,jk,jpzmi) + & 
     863                     trn(ji,jj,jk,jpzme) + trn(ji,jj,jk,jpdet) + trn(ji,jj,jk,jpdin) 
     864               loc_cycletot3(1) = loc_cycletot3(1) + ( zq1 * cvol(ji,jj,jk) ) 
     865               ! silicon 
     866               zq1 = trn(ji,jj,jk,jppds) + trn(ji,jj,jk,jpsil) 
     867               loc_cycletot3(2) = loc_cycletot3(2) + ( zq1 * cvol(ji,jj,jk) ) 
     868               ! iron 
     869               zq1 = ((trn(ji,jj,jk,jpphn) + trn(ji,jj,jk,jpphd) + trn(ji,jj,jk,jpzmi) + & 
     870                     trn(ji,jj,jk,jpzme) + trn(ji,jj,jk,jpdet)) * xrfn) + trn(ji,jj,jk,jpfer) 
     871               loc_cycletot3(3) = loc_cycletot3(3) + ( zq1 * cvol(ji,jj,jk) ) 
     872               ! carbon 
     873               zq1 = (trn(ji,jj,jk,jpphn) * xthetapn)  + (trn(ji,jj,jk,jpphd) * xthetapd)  +  & 
     874                     (trn(ji,jj,jk,jpzmi) * xthetazmi) + (trn(ji,jj,jk,jpzme) * xthetazme) +  & 
     875                     trn(ji,jj,jk,jpdtc) + trn(ji,jj,jk,jpdic) 
     876               loc_cycletot3(4) = loc_cycletot3(4) + ( zq1 * cvol(ji,jj,jk) ) 
     877               ! alkalinity 
     878               zq1 = trn(ji,jj,jk,jpalk) 
     879               loc_cycletot3(5) = loc_cycletot3(5) + ( zq1 * cvol(ji,jj,jk) ) 
     880               ! oxygen 
     881               zq1 = trn(ji,jj,jk,jpoxy) 
     882               loc_cycletot3(6) = loc_cycletot3(6) + ( zq1 * cvol(ji,jj,jk) ) 
     883            ENDDO 
     884         ENDDO 
     885      ENDDO 
     886      ! 
     887      ! sediment cells 
     888      loc_cycletot2(:) = 0._wp 
     889      loc_area         = 0._wp 
     890      jk = 1 
     891      DO jj = nldj,nlej 
     892         DO ji = nldi,nlei 
     893            loc_area = loc_area + (e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     894            ! nitrogen 
     895            zq1 = zn_sed_n(ji,jj) 
     896            loc_cycletot2(1) = loc_cycletot2(1) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     897            ! silicon 
     898            zq1 = zn_sed_si(ji,jj) 
     899            loc_cycletot2(2) = loc_cycletot2(2) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     900            ! iron 
     901            zq1 = zn_sed_fe(ji,jj) 
     902            loc_cycletot2(3) = loc_cycletot2(3) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     903            ! carbon 
     904            zq1 = zn_sed_c(ji,jj) + zn_sed_ca(ji,jj) 
     905            loc_cycletot2(4) = loc_cycletot2(4) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     906            ! alkalinity 
     907            zq1 = zn_sed_ca(ji,jj) * 2._wp 
     908            loc_cycletot2(5) = loc_cycletot2(5) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 
     909            ! skip oxygen 
     910         ENDDO 
     911      ENDDO 
     912     IF ( lwp ) WRITE(numout,*) '-- New Check --' 
     913      ! 
     914      !---------------------------------------------------------------------- 
     915      ! nitrogen 
     916      ! total tracer, and delta 
     917      zq1 = loc_cycletot3(1) 
     918      zq2 = loc_cycletot2(1) 
     919      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     920      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     921      zinvt      = zq1 + zq2 
     922      zdelta     = zinvt - cycletot2(1) 
     923      zratio     = 1.0e2 * zdelta / cycletot2(1) 
     924      ! 
     925      IF ( lwp ) WRITE(numout,9010) 'nitrogen', zq1, zq2, zinvt,   & 
     926         cycletot2(1), zdelta, zratio 
     927      IF ( lwp ) WRITE(numout,*) 
     928      ! 
     929      !---------------------------------------------------------------------- 
     930      ! silicon 
     931      zq1 = loc_cycletot3(2) 
     932      zq2 = loc_cycletot2(2) 
     933      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     934      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     935      zinvt      = zq1 + zq2 
     936      zdelta     = zinvt - cycletot2(2) 
     937      zratio     = 1.0e2 * zdelta / cycletot2(2) 
     938      ! 
     939      IF ( lwp ) WRITE(numout,9010) 'silicon', zq1, zq2, zinvt,   & 
     940         cycletot2(2), zdelta, zratio 
     941      IF ( lwp ) WRITE(numout,*) 
     942      ! 
     943      !---------------------------------------------------------------------- 
     944      ! iron 
     945      zq1 = loc_cycletot3(3) 
     946      zq2 = loc_cycletot2(3) 
     947      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     948      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     949      zinvt      = zq1 + zq2 
     950      zdelta     = zinvt - cycletot2(3) 
     951      zratio     = 1.0e2 * zdelta / cycletot2(3) 
     952      ! 
     953      IF ( lwp ) WRITE(numout,9010) 'iron', zq1, zq2, zinvt,   & 
     954         cycletot2(3), zdelta, zratio 
     955      IF ( lwp ) WRITE(numout,*) 
     956      ! 
     957      !---------------------------------------------------------------------- 
     958      ! carbon 
     959      zq1 = loc_cycletot3(4) 
     960      zq2 = loc_cycletot2(4) 
     961      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     962      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     963      zinvt      = zq1 + zq2 
     964      zdelta     = zinvt - cycletot2(4) 
     965      zratio     = 1.0e2 * zdelta / cycletot2(4) 
     966      ! 
     967      IF ( lwp ) WRITE(numout,9010) 'carbon', zq1, zq2, zinvt,   & 
     968         cycletot2(4), zdelta, zratio 
     969      IF ( lwp ) WRITE(numout,*) 
     970      ! 
     971      !---------------------------------------------------------------------- 
     972      ! alkalinity 
     973      zq1 = loc_cycletot3(5) 
     974      zq2 = loc_cycletot2(5) 
     975      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     976      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     977      zinvt      = zq1 + zq2 
     978      zdelta     = zinvt - cycletot2(5) 
     979      zratio     = 1.0e2 * zdelta / cycletot2(5) 
     980      ! 
     981      IF ( lwp ) WRITE(numout,9010) 'alkalinity', zq1, zq2, zinvt,   & 
     982         cycletot2(5), zdelta, zratio 
     983      IF ( lwp ) WRITE(numout,*) 
     984      ! 
     985      !---------------------------------------------------------------------- 
     986      ! oxygen 
     987      zq1 = loc_cycletot3(6) 
     988      zq2 = loc_cycletot2(6) 
     989      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     990      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     991      zinvt      = zq1 + zq2 
     992      zdelta     = zinvt - cycletot2(6) 
     993      zratio     = 1.0e2 * zdelta / cycletot2(6) 
     994      ! 
     995      IF ( lwp ) WRITE(numout,9010) 'oxygen', zq1, zq2, zinvt,   & 
     996         cycletot2(6), zdelta, zratio 
     997      IF ( lwp ) WRITE(numout,*) 
     998      ! 
     999      !--------------------------------------------------------------------- 
     1000      zq1 = loc_vol 
     1001      zq2 = loc_area 
     1002      IF( lk_mpp )   CALL mpp_sum( zq1 ) 
     1003      IF( lk_mpp )   CALL mpp_sum( zq2 ) 
     1004      IF ( lwp ) WRITE(numout,*) 
     1005      IF ( lwp ) WRITE(numout,*) ' check : cvol    : ', zq1 
     1006      IF ( lwp ) WRITE(numout,*) ' check : carea   : ', zq2 
     1007      IF ( lwp ) WRITE(numout,*) 
     1008      IF ( lwp ) CALL flush(numout) 
     1009 
     1010      !! 
     1011      !! 
     10129010  FORMAT(' element:',a10,                     & 
     1013             ' 3d sum:',e18.10,' 2d sum:',e18.10, & 
     1014             ' total:',e18.10,' initial:',e18.10, & 
     1015             ' delta:',e18.10,' %:',e18.10) 
     1016      ! 
     1017   END SUBROUTINE trc_rst_conserve  
     1018#endif 
     1019 
    6991020   SUBROUTINE trc_rst_tra_stat 
    7001021      !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.