Ignore:
Timestamp:
2018-02-16T18:31:39+01:00 (2 years ago)
Author:
dford
Message:

Add balancing code.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90

    r7803 r9331  
    2121   USE par_fabm                          ! FABM parameters 
    2222   USE par_trc,     ONLY: jptra          ! Tracer parameters 
     23   USE bioanalysis 
    2324 
    2425   IMPLICIT NONE 
     
    2728   PUBLIC asm_logchl_bal_ersem 
    2829 
     30   ! Default values for biological assimilation parameters 
     31   ! Should match Hemmings et al. (2008) 
     32   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor 
     33   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient 
     34   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate 
     35   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator 
     36   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator 
     37   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator 
     38   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator 
     39   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator 
     40   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator 
     41   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D 
     42   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction 
     43   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor 
     44   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor 
     45   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor 
     46   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor 
     47   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary) 
     48   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment 
     49   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation 
     50   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f. 
     51   ! 
     52   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics 
     53   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics 
     54   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing 
     55   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing 
     56   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD 
     57   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD 
     58   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD 
     59   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments 
     60 
    2961CONTAINS 
    3062 
    3163   SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, & 
    32       &                             k_maxchlinc, logchl_bkginc, logchl_balinc ) 
     64      &                             k_maxchlinc, logchl_bkginc, logchl_balinc,   & 
     65      &                             pgrow_avg_bkg, ploss_avg_bkg,           & 
     66      &                             phyt_avg_bkg, mld_max_bkg,              & 
     67      &                             ld_logchlbal, aincper ) 
    3368      !!--------------------------------------------------------------------------- 
    3469      !!                    ***  ROUTINE asm_logchl_bal_ersem  *** 
     
    5287      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,npfts)     :: logchl_bkginc 
    5388      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc 
     89      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth 
     90      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss 
     91      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto 
     92      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD 
     93      LOGICAL,  INTENT(in   )                               :: ld_logchlbal   ! Balancing y/n 
     94      REAL(wp), INTENT(in   )                               :: aincper 
    5495      !! 
    5596      INTEGER                      :: ji, jj, jk 
     
    5798      REAL(wp)                     :: chl_tot, chl_inc 
    5899      REAL(wp), DIMENSION(jpi,jpj) :: zmld 
     100      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters 
     101      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters 
     102      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state 
     103      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments 
     104      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics 
     105      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics 
     106      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices 
     107      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy 
     108      REAL(wp),                DIMENSION(jpi,jpj)           :: chlinctot 
     109      REAL(wp)                                              :: n2be_p         ! N:biomass for total phy 
     110      REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo 
     111      REAL(wp)                                              :: n2be_d         ! N:biomass for detritus 
     112      REAL(wp)                                              :: zfrac_n3n      ! Fraction of jp_fabm_n3n 
     113      REAL(wp)                                              :: zfrac_n4n      ! Fraction of jp_fabm_n4n 
     114      REAL(wp)                                              :: zfrac_r4n      ! Fraction of jp_fabm_r4n 
     115      REAL(wp)                                              :: zfrac_r6n      ! Fraction of jp_fabm_r6n 
     116      REAL(wp)                                              :: zfrac_r8n      ! Fraction of jp_fabm_r8n 
     117      REAL(wp)                                              :: zfrac_z4n      ! Fraction of jp_fabm_z4n 
     118      REAL(wp)                                              :: zfrac_z5n      ! Fraction of jp_fabm_z5n 
     119      REAL(wp)                                              :: zfrac_z6n      ! Fraction of jp_fabm_z6n 
     120      REAL(wp)                                              :: zfrac_p1n      ! Fraction of jp_fabm_p1n 
     121      REAL(wp)                                              :: zfrac_p2n      ! Fraction of jp_fabm_p2n 
     122      REAL(wp)                                              :: zfrac_p3n      ! Fraction of jp_fabm_p3n 
     123      REAL(wp)                                              :: zfrac_p4n      ! Fraction of jp_fabm_p4n 
     124      REAL(wp)                                              :: zrat_z4c_z4n 
     125      REAL(wp)                                              :: zrat_z5c_z5n 
     126      REAL(wp)                                              :: zrat_z5p_z5n 
     127      REAL(wp)                                              :: zrat_z6c_z6n 
     128      REAL(wp)                                              :: zrat_z6p_z6n 
     129      REAL(wp)                                              :: zrat_p1c_p1n 
     130      REAL(wp)                                              :: zrat_p1p_p1n 
     131      REAL(wp)                                              :: zrat_p1s_p1n 
     132      REAL(wp)                                              :: zrat_p2c_p2n 
     133      REAL(wp)                                              :: zrat_p2p_p2n 
     134      REAL(wp)                                              :: zrat_p3c_p3n 
     135      REAL(wp)                                              :: zrat_p3p_p3n 
     136      REAL(wp)                                              :: zrat_p4c_p4n 
     137      REAL(wp)                                              :: zrat_p4p_p4n 
    59138      !!--------------------------------------------------------------------------- 
    60139 
     
    176255         zmld(:,:) = hmlpt(:,:) 
    177256      END SELECT 
     257 
     258      IF ( ld_logchlbal ) THEN 
     259       
     260         ! Set up model parameters to be passed into Hemmings balancing routine. 
     261         ! For now these are hardwired to the standard HadOCC parameter values 
     262         ! as this is what the scheme was developed for. 
     263         ! Obviously, HadOCC and ERSEM are rather different models, so this 
     264         ! isn't ideal, but there's not always direct analogues between the two 
     265         ! parameter sets, so it's the easiest way to get something running. 
     266         ! In the longer term, some serious MarMOT-based development is required. 
     267         modparm(1)  = 0.1                               ! grow_sat 
     268         modparm(2)  = 2.0                               ! psmax 
     269         modparm(3)  = 0.845                             ! par 
     270         modparm(4)  = 0.02                              ! alpha 
     271         modparm(5)  = 0.05                              ! resp_rate 
     272         modparm(6)  = 0.05                              ! pmort_rate 
     273         modparm(7)  = 0.01                              ! phyto_min 
     274         modparm(8)  = 0.05                              ! z_mort_1 
     275         modparm(9)  = 1.0                               ! z_mort_2 
     276         modparm(10) = 6.625                             ! c2n_p 
     277         modparm(11) = 5.625                             ! c2n_z 
     278         modparm(12) = 7.5                               ! c2n_d 
     279         modparm(13) = 0.01                              ! graze_threshold 
     280         modparm(14) = 2.0                               ! holling_coef 
     281         modparm(15) = 0.5                               ! graze_sat 
     282         modparm(16) = 2.0                               ! graze_max 
     283 
     284         ! Set up assimilation parameters to be passed into balancing routine 
     285         ! Not sure what assimparm(1) is meant to be, but it doesn't get used 
     286         assimparm(2)  = balnutext 
     287         assimparm(3)  = balnutmin 
     288         assimparm(4)  = r 
     289         assimparm(5)  = beta_g 
     290         assimparm(6)  = beta_l 
     291         assimparm(7)  = beta_m 
     292         assimparm(8)  = a_g 
     293         assimparm(9)  = a_l 
     294         assimparm(10) = a_m 
     295         assimparm(11) = zfracb0 
     296         assimparm(12) = zfracb1 
     297         assimparm(13) = qrfmax 
     298         assimparm(14) = qafmax 
     299         assimparm(15) = zrfmax 
     300         assimparm(16) = zafmax 
     301         assimparm(17) = prfmax 
     302         assimparm(18) = incphymin 
     303         assimparm(19) = integnstep 
     304         assimparm(20) = pthreshold 
     305 
     306         ! Set up external tracer indices array bstate 
     307         i_tracer(1) = 1   ! nutrient 
     308         i_tracer(2) = 2   ! phytoplankton 
     309         i_tracer(3) = 3   ! zooplankton 
     310         i_tracer(4) = 4   ! detritus 
     311         i_tracer(5) = 5   ! DIC 
     312         i_tracer(6) = 6   ! Alkalinity 
     313 
     314         ! Set background state 
     315         bstate(:,:,:,i_tracer(1)) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) + trn(:,:,:,jp_fabm_m1+jp_fabm_n4n) 
     316         bstate(:,:,:,i_tracer(2)) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1n) + trn(:,:,:,jp_fabm_m1+jp_fabm_p2n) + & 
     317            &                        trn(:,:,:,jp_fabm_m1+jp_fabm_p3n) + trn(:,:,:,jp_fabm_m1+jp_fabm_p4n) 
     318         ! Z4c needs converting by qnc, hardwire for now 
     319         bstate(:,:,:,i_tracer(3)) = (trn(:,:,:,jp_fabm_m1+jp_fabm_z4c) * 0.0126 ) + & 
     320            &                        trn(:,:,:,jp_fabm_m1+jp_fabm_z5n) + trn(:,:,:,jp_fabm_m1+jp_fabm_z6n) 
     321         bstate(:,:,:,i_tracer(4)) = trn(:,:,:,jp_fabm_m1+jp_fabm_r4n) + trn(:,:,:,jp_fabm_m1+jp_fabm_r6n) + & 
     322            &                        trn(:,:,:,jp_fabm_m1+jp_fabm_r8n) 
     323         bstate(:,:,:,i_tracer(5)) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) 
     324         bstate(:,:,:,i_tracer(6)) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3a) 
     325 
     326         ! Calculate carbon to chlorophyll ratio for combined phytoplankton 
     327         ! and nitrogen to biomass equivalent for PZD 
     328         ! Need a single number, so base on HadOCC 
     329         cchl_p(:,:) = 0.0 
     330         DO jj = 1, jpj 
     331            DO ji = 1, jpi 
     332               IF ( ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + & 
     333                  &   trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) ) .GT. 0.0 ) THEN 
     334                  cchl_p(ji,jj) = ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_p1c)  + trn(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) + & 
     335                     &              trn(ji,jj,1,jp_fabm_m1+jp_fabm_p3c)  + trn(ji,jj,1,jp_fabm_m1+jp_fabm_p4c) ) /  & 
     336                     &            ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + & 
     337                     &              trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) ) 
     338               ENDIF 
     339            END DO 
     340         END DO 
     341         n2be_p = ( 14.01 + ( 12.01 * 6.625 ) ) ! / ( 14.01 + ( 12.01 * 6.625 ) ) 
     342         n2be_z = ( 14.01 + ( 12.01 * 5.625 ) ) ! / ( 14.01 + ( 12.01 * 6.625 ) ) 
     343         n2be_d = ( 14.01 + ( 12.01 * 7.5 )   ) ! / ( 14.01 + ( 12.01 * 6.625 ) ) 
     344 
     345         chlinctot(:,:) = logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl1) + logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
     346            &             logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl3) + logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     347 
     348         ! Call nitrogen balancing routine 
     349         CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   & 
     350            &               n2be_p, n2be_z, n2be_d, assimparm,                      & 
     351            &               INT(aincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       & 
     352            &               zmld(:,:), mld_max_bkg(:,:), chlinctot(:,:), cchl_p(:,:), & 
     353            &               nbal_active, phyt_avg_bkg(:,:),                         & 
     354            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      & 
     355            &               subsurf_active, deepneg_active,                         & 
     356            &               deeppos_active, nutprof_active,                         & 
     357            &               bstate, outincs,                                        & 
     358            &               diag_active, diag,                                      & 
     359            &               diag_fulldepth_active, diag_fulldepth ) 
     360          
     361         ! Loop over each grid point partioning the increments 
     362         DO jk = 1, jpk 
     363            DO jj = 1, jpj 
     364               DO ji = 1, jpi 
     365 
     366                  ! Nitrogen phytoplankton from balancing scheme 
     367                  ! Split according to current ratios [ChlTot] or assimilation [PFTs] 
     368                  ! Update carbon, phosphorus and silicon according to current ratios 
     369                  ! Already have chlorophyll 
     370                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) > 0.0 ) .AND. & 
     371                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) > 0.0 ) .AND. & 
     372                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) > 0.0 ) .AND. & 
     373                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) > 0.0 ) ) THEN 
     374 
     375                     IF ( ld_logchlpftinc ) THEN 
     376                        IF ( (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + & 
     377                           &  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + & 
     378                           &  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + & 
     379                           &  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)) > 0.0 ) THEN 
     380                           zfrac_p1n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / & 
     381                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + & 
     382                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + & 
     383                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + & 
     384                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)) 
     385                           zfrac_p2n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / & 
     386                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + & 
     387                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + & 
     388                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + & 
     389                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)) 
     390                           zfrac_p3n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / & 
     391                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + & 
     392                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + & 
     393                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + & 
     394                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)) 
     395                           zfrac_p4n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / & 
     396                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + & 
     397                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + & 
     398                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + & 
     399                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)) 
     400                        ELSE 
     401                           zfrac_p1n = 0.25 
     402                           zfrac_p2n = 0.25 
     403                           zfrac_p3n = 0.25 
     404                           zfrac_p4n = 0.25 
     405                        ENDIF 
     406                     ELSE 
     407                        zfrac_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) / & 
     408                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + & 
     409                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + & 
     410                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + & 
     411                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)) 
     412                        zfrac_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) / & 
     413                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + & 
     414                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + & 
     415                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + & 
     416                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)) 
     417                        zfrac_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) / & 
     418                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + & 
     419                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + & 
     420                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + & 
     421                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)) 
     422                        zfrac_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) / & 
     423                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + & 
     424                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + & 
     425                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + & 
     426                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)) 
     427                     ENDIF 
     428                      
     429                     zrat_p1c_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) 
     430                     zrat_p1p_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) 
     431                     zrat_p1s_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) 
     432                     zrat_p2c_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) 
     433                     zrat_p2p_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) 
     434                     zrat_p3c_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) 
     435                     zrat_p3p_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) 
     436                     zrat_p4c_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) 
     437                     zrat_p4p_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) 
     438                      
     439                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p1n 
     440                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p2n 
     441                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p3n 
     442                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p4n 
     443                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1c_p1n 
     444                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1p_p1n 
     445                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1s_p1n 
     446                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) * zrat_p2c_p2n 
     447                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) * zrat_p2p_p2n 
     448                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) * zrat_p3c_p3n 
     449                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) * zrat_p3p_p3n 
     450                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) * zrat_p4c_p4n 
     451                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) * zrat_p4p_p4n 
     452 
     453                  ENDIF 
     454 
     455                  ! Nitrogen nutrient from balancing scheme 
     456                  ! Split between nitrate and ammonium according to current ratios 
     457                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) > 0.0 ) .AND. & 
     458                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) > 0.0 ) ) THEN 
     459                     zfrac_n3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) / & 
     460                        &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) + trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n)) 
     461                     zfrac_n4n = 1.0 - zfrac_n3n 
     462                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n3n 
     463                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n4n 
     464                  ENDIF 
     465 
     466                  ! Nitrogen zooplankton from balancing scheme 
     467                  ! Split according to current ratios 
     468                  ! Update carbon and phosphorus according to current ratios 
     469                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) > 0.0 ) .AND. & 
     470                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) > 0.0 ) .AND. & 
     471                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) > 0.0 ) ) THEN 
     472                     zfrac_z4n = (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) / & 
     473                        &        ((trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) + & 
     474                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + & 
     475                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n)) 
     476                     zfrac_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) / & 
     477                        &        ((trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) + & 
     478                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + & 
     479                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n)) 
     480                     zfrac_z6n = 1.0 - zfrac_z4n - zfrac_z5n 
     481                     zrat_z4c_z4n = 1.0 / 0.0126 
     482                     zrat_z5c_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) 
     483                     zrat_z5p_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) 
     484                     zrat_z6c_z6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) 
     485                     zrat_z6p_z6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) 
     486                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z5n 
     487                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z6n 
     488                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z4n * zrat_z4c_z4n 
     489                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) * zrat_z5c_z5n 
     490                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) * zrat_z6c_z6n 
     491                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) * zrat_z5p_z5n 
     492                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) * zrat_z6p_z6n 
     493                  ENDIF 
     494 
     495                  ! Nitrogen detritus from balancing scheme 
     496                  ! Split according to current ratios 
     497                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) > 0.0 ) .AND. & 
     498                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) > 0.0 ) .AND. & 
     499                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) > 0.0 ) ) THEN 
     500                     zfrac_r4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) / & 
     501                        &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + & 
     502                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + & 
     503                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n)) 
     504                     zfrac_r6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) / & 
     505                        &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + & 
     506                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + & 
     507                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n)) 
     508                     zfrac_r8n = 1.0 - zfrac_r4n - zfrac_r6n 
     509                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r4n 
     510                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r6n 
     511                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r8n 
     512                  ENDIF 
     513 
     514                  ! DIC straight from balancing scheme 
     515                  logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3c) = outincs(ji,jj,jk,i_tracer(5)) 
     516 
     517                  ! Alkalinity straight from balancing scheme 
     518                  logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3a) = outincs(ji,jj,jk,i_tracer(6)) 
     519       
     520               END DO 
     521            END DO 
     522         END DO 
     523       
     524      ENDIF 
    178525      ! 
    179       ! Now set MLD to bottom of a level and propagate incs equally through mixed layer 
     526      ! Now set MLD to bottom of a level and propagate chlorophyll incs equally through mixed layer 
     527      ! If balancing, should really relate this back to phytoplankton, but stick with this for now 
    180528      ! 
    181529      DO jj = 1, jpj 
     
    201549      END DO 
    202550 
    203       ! Multivariate balancing forthcoming... 
    204  
    205551   END SUBROUTINE asm_logchl_bal_ersem 
    206552 
Note: See TracChangeset for help on using the changeset viewer.