Ignore:
Timestamp:
2015-07-15T17:46:12+02:00 (5 years ago)
Author:
andrewryan
Message:

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90

    r5034 r5600  
    2424   USE p4zsed          !  Sedimentation 
    2525   USE p4zint          !  time interpolation 
     26   USE p4zrem          !  remineralisation 
    2627   USE iom             !  I/O manager 
    2728   USE trd_oce         !  Ocean trends variables 
     
    3637   PUBLIC   p4z_sms        ! called in p4zsms.F90 
    3738 
    38    REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget 
    39    INTEGER ::  numco2, numnut  !: logical unit for co2 budget 
    40  
     39   REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget, po4budget 
     40   REAL(wp) :: xfact1, xfact2, xfact3 
     41   INTEGER ::  numco2, numnut, numnit  !: logical unit for co2 budget 
     42 
     43   !!* Array used to indicate negative tracer values 
     44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr     !: ??? 
     45 
     46 
     47   !! * Substitutions 
     48#  include "top_substitute.h90" 
    4149   !!---------------------------------------------------------------------- 
    4250   !! NEMO/TOP 3.3 , NEMO Consortium (2010) 
     
    6169      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index       
    6270      !! 
    63       INTEGER ::   jnt, jn, jl 
     71      INTEGER ::   ji, jj, jk, jnt, jn, jl 
     72      REAL(wp) ::  ztra 
     73#if defined key_kriest 
     74      REAL(wp) ::  zcoef1, zcoef2 
     75#endif 
    6476      CHARACTER (len=25) :: charout 
    65       REAL(wp), POINTER, DIMENSION(:,:,:,:)  :: ztrdpis 
    6677      !!--------------------------------------------------------------------- 
    6778      ! 
    6879      IF( nn_timing == 1 )  CALL timing_start('p4z_sms') 
    6980      ! 
    70       IF( l_trdtrc )  THEN 
    71          CALL wrk_alloc( jpi, jpj, jpk, jp_pisces, ztrdpis )  
    72          DO jn = 1, jp_pisces 
    73             jl = jn + jp_pcs0 - 1 
    74             ztrdpis(:,:,:,jn) = trn(:,:,:,jl) 
    75          ENDDO 
    76       ENDIF 
    77       ! 
    7881      IF( kt == nittrc000 ) THEN 
     82        ! 
     83        ALLOCATE( xnegtr(jpi,jpj,jpk) ) 
    7984        ! 
    8085        CALL p4z_che                              ! initialize the chemical constants 
     
    8893      IF( ln_pisdmp .AND. MOD( kt - nn_dttrc, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt )      ! Relaxation of some tracers 
    8994      ! 
     95      !                                                                    !   set time step size (Euler/Leapfrog) 
     96      IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN   ;    rfact = rdttrc(1)     !  at nittrc000 
     97      ELSEIF( kt <= nittrc000 + nn_dttrc )                          THEN   ;    rfact = 2. * rdttrc(1)   ! at nittrc000 or nittrc000+nn_dttrc (Leapfrog) 
     98      ENDIF 
     99      ! 
     100      IF( ( ln_top_euler .AND. kt == nittrc000 )  .OR. ( .NOT.ln_top_euler .AND. kt <= nittrc000 + nn_dttrc ) ) THEN 
     101         rfactr  = 1. / rfact 
     102         rfact2  = rfact / FLOAT( nrdttrc ) 
     103         rfact2r = 1. / rfact2 
     104         xstep = rfact2 / rday         ! Time step duration for biology 
     105         IF(lwp) WRITE(numout,*)  
     106         IF(lwp) WRITE(numout,*) '    Passive Tracer  time step    rfact  = ', rfact, ' rdt = ', rdttra(1) 
     107         IF(lwp) write(numout,*) '    PISCES  Biology time step    rfact2 = ', rfact2 
     108         IF(lwp) WRITE(numout,*) 
     109      ENDIF 
     110 
     111      IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN 
     112         DO jn = jp_pcs0, jp_pcs1              !   SMS on tracer without Asselin time-filter 
     113            trb(:,:,:,jn) = trn(:,:,:,jn) 
     114         END DO 
     115      ENDIF 
     116      ! 
    90117      IF( ndayflxtr /= nday_year ) THEN      ! New days 
    91118         ! 
     
    105132      DO jnt = 1, nrdttrc          ! Potential time splitting if requested 
    106133         ! 
    107          CALL p4z_bio (kt, jnt)    ! Biology 
    108          CALL p4z_sed (kt, jnt)    ! Sedimentation 
    109          ! 
     134         CALL p4z_bio( kt, jnt )   ! Biology 
     135         CALL p4z_sed( kt, jnt )   ! Sedimentation 
     136         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation 
     137         CALL p4z_flx( kt, jnt )   ! Compute surface fluxes 
     138         ! 
     139         xnegtr(:,:,:) = 1.e0 
    110140         DO jn = jp_pcs0, jp_pcs1 
    111             trb(:,:,:,jn) = trn(:,:,:,jn) 
    112          ENDDO 
    113          ! 
     141            DO jk = 1, jpk 
     142               DO jj = 1, jpj 
     143                  DO ji = 1, jpi 
     144                     IF( ( trb(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) < 0.e0 ) THEN 
     145                        ztra             = ABS( trb(ji,jj,jk,jn) ) / ( ABS( tra(ji,jj,jk,jn) ) + rtrn ) 
     146                        xnegtr(ji,jj,jk) = MIN( xnegtr(ji,jj,jk),  ztra ) 
     147                     ENDIF 
     148                 END DO 
     149               END DO 
     150            END DO 
     151         END DO 
     152         !                                ! where at least 1 tracer concentration becomes negative 
     153         !                                !  
     154         DO jn = jp_pcs0, jp_pcs1 
     155           trb(:,:,:,jn) = trb(:,:,:,jn) + xnegtr(:,:,:) * tra(:,:,:,jn) 
     156         END DO 
     157        ! 
     158         DO jn = jp_pcs0, jp_pcs1 
     159            tra(:,:,:,jn) = 0._wp 
     160         END DO 
     161         ! 
     162         IF( ln_top_euler ) THEN 
     163            DO jn = jp_pcs0, jp_pcs1 
     164               trn(:,:,:,jn) = trb(:,:,:,jn) 
     165            END DO 
     166         ENDIF 
    114167      END DO 
    115168 
    116       IF( l_trdtrc )  THEN 
    117          DO jn = 1, jp_pisces 
    118             jl = jn + jp_pcs0 - 1 
    119             ztrdpis(:,:,:,jn) = ( ztrdpis(:,:,:,jn) - trn(:,:,:,jl) ) * rfact2r 
    120          ENDDO 
    121       ENDIF 
    122       CALL p4z_lys( kt )             ! Compute CaCO3 saturation 
    123       CALL p4z_flx( kt )             ! Compute surface fluxes 
    124  
    125       DO jn = jp_pcs0, jp_pcs1 
    126         CALL lbc_lnk( trn(:,:,:,jn), 'T', 1. ) 
    127         CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. ) 
    128         CALL lbc_lnk( tra(:,:,:,jn), 'T', 1. ) 
     169#if defined key_kriest 
     170      !  
     171      zcoef1 = 1.e0 / xkr_massp  
     172      zcoef2 = 1.e0 / xkr_massp / 1.1 
     173      DO jk = 1,jpkm1 
     174         trb(:,:,jk,jpnum) = MAX(  trb(:,:,jk,jpnum), trb(:,:,jk,jppoc) * zcoef1 / xnumm(jk)  ) 
     175         trb(:,:,jk,jpnum) = MIN(  trb(:,:,jk,jpnum), trb(:,:,jk,jppoc) * zcoef2              ) 
    129176      END DO 
    130177      ! 
     178#endif 
     179      ! 
     180      ! 
     181      IF( l_trdtrc ) THEN 
     182         DO jn = jp_pcs0, jp_pcs1 
     183           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends 
     184         END DO 
     185      END IF 
     186      ! 
    131187      IF( lk_sed ) THEN  
    132188         ! 
     
    134190         ! 
    135191         DO jn = jp_pcs0, jp_pcs1 
    136            CALL lbc_lnk( trn(:,:,:,jn), 'T', 1. ) 
     192           CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. ) 
    137193         END DO 
    138194         ! 
     
    141197      IF( lrst_trc )  CALL p4z_rst( kt, 'WRITE' )  !* Write PISCES informations in restart file  
    142198      ! 
    143       IF( l_trdtrc ) THEN 
    144          DO jn = 1, jp_pisces 
    145             jl = jn + jp_pcs0 - 1 
    146              ztrdpis(:,:,:,jn) = ztrdpis(:,:,:,jn) + tra(:,:,:,jl) 
    147              CALL trd_trc( ztrdpis(:,:,:,jn), jn, jptra_sms, kt )   ! save trends 
    148           END DO 
    149           CALL wrk_dealloc( jpi, jpj, jpk, jp_pisces, ztrdpis )  
    150       END IF 
    151       ! 
     199 
    152200      IF( lk_iomput .OR. ln_check_mass )  CALL p4z_chk_mass( kt ) ! Mass conservation checking 
    153201 
     
    280328               ztmas   = tmask(ji,jj,jk) 
    281329               ztmas1  = 1. - tmask(ji,jj,jk) 
    282                zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  ) 
    283                zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1 
    284                zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk ) 
     330               zcaralk = trb(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  ) 
     331               zco3    = ( zcaralk - trb(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1 
     332               zbicarb = ( 2. * trb(ji,jj,jk,jpdic) - zcaralk ) 
    285333               hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1 
    286334            END DO 
     
    361409      REAL(wp) ::  silmean = 91.51     ! mean value of silicate 
    362410      ! 
    363       REAL(wp) :: zarea, zalksum, zpo4sum, zno3sum, zsilsum 
     411      REAL(wp) :: zarea, zalksumn, zpo4sumn, zno3sumn, zsilsumn 
     412      REAL(wp) :: zalksumb, zpo4sumb, zno3sumb, zsilsumb 
    364413      !!--------------------------------------------------------------------- 
    365414 
     
    374423         zarea          = 1._wp / glob_sum( cvol(:,:,:) ) * 1e6               
    375424 
    376          zalksum = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea 
    377          zpo4sum = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r 
    378          zno3sum = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3 
    379          zsilsum = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea 
     425         zalksumn = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea 
     426         zpo4sumn = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r 
     427         zno3sumn = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3 
     428         zsilsumn = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea 
    380429  
    381          IF(lwp) WRITE(numout,*) '       TALK mean : ', zalksum 
    382          trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksum 
    383  
    384          IF(lwp) WRITE(numout,*) '       PO4  mean : ', zpo4sum 
    385          trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sum 
    386  
    387          IF(lwp) WRITE(numout,*) '       NO3  mean : ', zno3sum 
    388          trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sum 
    389  
    390          IF(lwp) WRITE(numout,*) '       SiO3 mean : ', zsilsum 
    391          trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsum ) 
    392          ! 
    393       ENDIF 
    394  
     430         IF(lwp) WRITE(numout,*) '       TALKN mean : ', zalksumn 
     431         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksumn 
     432 
     433         IF(lwp) WRITE(numout,*) '       PO4N  mean : ', zpo4sumn 
     434         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sumn 
     435 
     436         IF(lwp) WRITE(numout,*) '       NO3N  mean : ', zno3sumn 
     437         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sumn 
     438 
     439         IF(lwp) WRITE(numout,*) '       SiO3N mean : ', zsilsumn 
     440         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsumn ) 
     441         ! 
     442         ! 
     443         IF( .NOT. ln_top_euler ) THEN 
     444            zalksumb = glob_sum( trb(:,:,:,jptal) * cvol(:,:,:)  ) * zarea 
     445            zpo4sumb = glob_sum( trb(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r 
     446            zno3sumb = glob_sum( trb(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3 
     447            zsilsumb = glob_sum( trb(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea 
     448  
     449            IF(lwp) WRITE(numout,*) ' ' 
     450            IF(lwp) WRITE(numout,*) '       TALKB mean : ', zalksumb 
     451            trb(:,:,:,jptal) = trb(:,:,:,jptal) * alkmean / zalksumb 
     452 
     453            IF(lwp) WRITE(numout,*) '       PO4B  mean : ', zpo4sumb 
     454            trb(:,:,:,jppo4) = trb(:,:,:,jppo4) * po4mean / zpo4sumb 
     455 
     456            IF(lwp) WRITE(numout,*) '       NO3B  mean : ', zno3sumb 
     457            trb(:,:,:,jpno3) = trb(:,:,:,jpno3) * no3mean / zno3sumb 
     458 
     459            IF(lwp) WRITE(numout,*) '       SiO3B mean : ', zsilsumb 
     460            trb(:,:,:,jpsil) = MIN( 400.e-6,trb(:,:,:,jpsil) * silmean / zsilsumb ) 
     461        ENDIF 
     462        ! 
     463      ENDIF 
     464        ! 
    395465   END SUBROUTINE p4z_dmp 
    396466 
     
    404474      !!--------------------------------------------------------------------- 
    405475      ! 
    406       INTEGER , INTENT( in ) ::   kt      ! ocean time-step index       
    407       REAL(wp)               ::  zfact        
    408       !! 
     476      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index       
     477      REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot 
     478      CHARACTER(LEN=100)   ::   cltxt 
     479      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol 
     480      INTEGER :: jk 
     481      !!---------------------------------------------------------------------- 
     482 
     483      ! 
    409484      !!--------------------------------------------------------------------- 
    410485 
     
    413488            CALL ctl_opn( numco2, 'carbon.budget'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea ) 
    414489            CALL ctl_opn( numnut, 'nutrient.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea ) 
     490            CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea ) 
     491            xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr 
     492            xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr 
     493            xfact3 = 1.e+3 * rfact2r * rno3   ! conversion molC/l/kt ----> molN/m3/s 
     494            cltxt='time-step   Alkalinity        Nitrate        Phosphorus         Silicate           Iron' 
     495            IF( lwp ) WRITE(numnut,*)  TRIM(cltxt) 
     496            IF( lwp ) WRITE(numnut,*)  
    415497         ENDIF 
    416498      ENDIF 
    417499 
     500      ! 
    418501      IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN 
    419502         !   Compute the budget of NO3, ALK, Si, Fer 
     
    431514      ENDIF 
    432515      ! 
    433       IF( iom_use( "psiltot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN 
     516      IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN 
     517         po4budget = glob_sum( (   trn(:,:,:,jppo4)                     & 
     518            &                    + trn(:,:,:,jpphy) + trn(:,:,:,jpdia)  & 
     519            &                    + trn(:,:,:,jpzoo) + trn(:,:,:,jpmes)  & 
     520            &                    + trn(:,:,:,jppoc)                     & 
     521#if ! defined key_kriest 
     522            &                    + trn(:,:,:,jpgoc)                     & 
     523#endif 
     524            &                    + trn(:,:,:,jpdoc)                     ) * cvol(:,:,:)  ) 
     525         po4budget = po4budget / areatot 
     526         CALL iom_put( "ppo4tot", po4budget ) 
     527      ENDIF 
     528      ! 
     529      IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN 
    434530         silbudget = glob_sum( (   trn(:,:,:,jpsil) + trn(:,:,:,jpgsi)  & 
    435531            &                    + trn(:,:,:,jpdsi)                     ) * cvol(:,:,:)  ) 
     
    439535      ENDIF 
    440536      ! 
    441       IF( iom_use( "palktot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN 
     537      IF( iom_use( "palktot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN 
    442538         alkbudget = glob_sum( (   trn(:,:,:,jpno3) * rno3              & 
    443539            &                    + trn(:,:,:,jptal)                     & 
     
    448544      ENDIF 
    449545      ! 
    450       IF( iom_use( "pfertot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN 
     546      IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN 
    451547         ferbudget = glob_sum( (   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe)  & 
    452548            &                    + trn(:,:,:,jpdfe)                     & 
     
    462558      ENDIF 
    463559      ! 
     560 
     561      ! Global budget of N SMS : denitrification in the water column and in the sediment 
     562      !                          nitrogen fixation by the diazotrophs 
     563      ! -------------------------------------------------------------------------------- 
     564      IF( iom_use( "tnfix" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN 
     565         znitrpottot  = glob_sum ( nitrpot(:,:,:) * nitrfix * cvol(:,:,:) ) 
     566         CALL iom_put( "tnfix"  , znitrpottot * 1.e+3 * rno3 )  ! Global  nitrogen fixation molC/l  to molN/m3  
     567      ENDIF 
     568      ! 
     569      IF( iom_use( "tdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN 
     570         zrdenittot   = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) ) 
     571         CALL iom_put( "tdenit"  , zrdenittot * 1.e+3 * rno3 )  ! Total denitrification molC/l to molN/m3  
     572      ENDIF 
     573      ! 
     574      IF( iom_use( "Sdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN 
     575         zsdenittot   = glob_sum ( sdenit(:,:) * e1e2t(:,:) ) 
     576         CALL iom_put( "Sdenit", sdenit(:,:) * xfact3 * tmask(:,:,1) )  ! Nitrate reduction in the sediments 
     577      ENDIF 
     578 
    464579      IF( ln_check_mass .AND. kt == nitend ) THEN   ! Compute the budget of NO3, ALK, Si, Fer 
    465          zfact = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/year 
    466580         t_atm_co2_flx  = t_atm_co2_flx / glob_sum( e1e2t(:,:) ) 
    467          t_oce_co2_flx  = t_oce_co2_flx         * zfact * (-1 ) 
    468          tpp            = tpp           * 1000. * zfact 
    469          t_oce_co2_exp  = t_oce_co2_exp * 1000. * zfact 
     581         t_oce_co2_flx  = t_oce_co2_flx         * xfact1 * (-1 ) 
     582         tpp            = tpp           * 1000. * xfact1 
     583         t_oce_co2_exp  = t_oce_co2_exp * 1000. * xfact1 
    470584         IF( lwp ) WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp 
    471          IF( lwp ) WRITE(numnut,9500) ndastp, alkbudget        * 1.e+06, & 
     585         IF( lwp ) WRITE(numnut,9100) ndastp, alkbudget        * 1.e+06, & 
    472586             &                                no3budget * rno3 * 1.e+06, & 
     587             &                                po4budget * po4r * 1.e+06, & 
    473588             &                                silbudget        * 1.e+06, & 
    474589             &                                ferbudget        * 1.e+09 
     590         ! 
     591         IF( lwp ) WRITE(numnit,9200) ndastp, znitrpottot * xfact2  , & 
     592         &                             zrdenittot  * xfact2  , & 
     593         &                             zsdenittot  * xfact2 
     594 
    475595      ENDIF 
    476596      ! 
    477597 9000  FORMAT(i8,f10.5,e18.10,f10.5,f10.5) 
    478  9500  FORMAT(i8,4e18.10) 
     598 9100  FORMAT(i8,5e18.10) 
     599 9200  FORMAT(i8,3f10.5) 
     600 
    479601       ! 
    480602   END SUBROUTINE p4z_chk_mass 
Note: See TracChangeset for help on using the changeset viewer.