New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11990 for branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto3dbal_medusa.F90 – NEMO

Ignore:
Timestamp:
2019-11-27T16:43:05+01:00 (4 years ago)
Author:
dford
Message:

Get the nitrogen balancing working with 3D chlorophyll increments.

File:
1 copied

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto3dbal_medusa.F90

    r11958 r11990  
    1 MODULE asmphyto2dbal_medusa 
     1MODULE asmphyto3dbal_medusa 
    22   !!====================================================================== 
    33   !!                       ***  MODULE asmphyto2dbal_medusa  *** 
     
    3434   PRIVATE                    
    3535 
    36    PUBLIC asm_phyto2d_bal_medusa 
     36   PUBLIC asm_phyto3d_bal_medusa 
    3737 
    3838   ! Default values for biological assimilation parameters 
     
    6969CONTAINS 
    7070 
    71    SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      & 
    72       &                               pinc_chltot,                    & 
    73       &                               ld_chldia,                      & 
    74       &                               pinc_chldia,                    & 
    75       &                               ld_chlnon,                      & 
    76       &                               pinc_chlnon,                    & 
    77       &                               ld_phytot,                      & 
    78       &                               pinc_phytot,                    & 
    79       &                               ld_phydia,                      & 
    80       &                               pinc_phydia,                    & 
    81       &                               ld_phynon,                      & 
    82       &                               pinc_phynon,                    & 
     71   SUBROUTINE asm_phyto3d_bal_medusa( pinc_chltot,                    & 
    8372      &                               pincper,                        & 
    84       &                               p_maxchlinc, ld_phytobal, pmld, & 
     73      &                               p_maxchlinc,                   & 
    8574      &                               pgrow_avg_bkg, ploss_avg_bkg,   & 
    86       &                               phyt_avg_bkg, mld_max_bkg,      & 
    87       &                               tracer_bkg, phyto2d_balinc ) 
     75      &                               phyt_avg_bkg,                   & 
     76      &                               tracer_bkg, phyto3d_balinc ) 
    8877      !!--------------------------------------------------------------------------- 
    8978      !!                    ***  ROUTINE asm_phyto2d_bal_medusa  *** 
     
    10190      !!--------------------------------------------------------------------------- 
    10291      !! 
    103       LOGICAL,  INTENT(in   )                               :: ld_chltot      ! Assim chltot y/n 
    104       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chltot    ! chltot increments 
    105       LOGICAL,  INTENT(in   )                               :: ld_chldia      ! Assim chldia y/n 
    106       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chldia    ! chldia increments 
    107       LOGICAL,  INTENT(in   )                               :: ld_chlnon      ! Assim chlnon y/n 
    108       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chlnon    ! chlnon increments 
    109       LOGICAL,  INTENT(in   )                               :: ld_phytot      ! Assim phytot y/n 
    110       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phytot    ! phytot increments 
    111       LOGICAL,  INTENT(in   )                               :: ld_phydia      ! Assim phydia y/n 
    112       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phydia    ! phydia increments 
    113       LOGICAL,  INTENT(in   )                               :: ld_phynon      ! Assim phynon y/n 
    114       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phynon    ! phynon increments 
     92      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk)       :: pinc_chltot    ! chltot increments 
    11593      REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period 
    11694      REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment 
    117       LOGICAL,  INTENT(in   )                               :: ld_phytobal    ! Balancing y/n 
    118       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld           ! Mixed layer depth 
    119       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth 
    120       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss 
    121       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto 
    122       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD 
     95      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)       :: pgrow_avg_bkg  ! Avg phyto growth 
     96      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)       :: ploss_avg_bkg  ! Avg phyto loss 
     97      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)       :: phyt_avg_bkg   ! Avg phyto 
    12398      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables 
    124       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments 
     99      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto3d_balinc ! Balancing increments 
    125100      !! 
    126101      INTEGER                                               :: ji, jj, jk, jn ! Loop counters 
     
    144119      REAL(wp)                                              :: zrat_pds_chd   ! Ratio of jppds:jpchd 
    145120      REAL(wp)                                              :: zrat_dtc_det   ! Ratio of jpdtc:jpdet 
    146       REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy 
     121      REAL(wp),                DIMENSION(jpi,jpj,jpk)       :: cchl_p         ! C:Chl for total phy 
     122      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p_2d      ! C:Chl for total phy 
    147123      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters 
    148124      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters 
    149125      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state 
     126      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: bstate_2d      ! Background state 
    150127      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments 
     128      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: outincs_2d     ! Balancing increments 
    151129      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics 
    152       REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics 
     130      REAL(wp),                DIMENSION(jpi,jpj,1,22)      :: diag_fulldepth ! Full-depth diagnostics 
     131      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chltot_2d 
     132      REAL(wp),                DIMENSION(jpi,jpj)           :: pgrow_avg_bkg_2d 
     133      REAL(wp),                DIMENSION(jpi,jpj)           :: ploss_avg_bkg_2d 
     134      REAL(wp),                DIMENSION(jpi,jpj)           :: phyt_avg_bkg_2d 
     135      REAL(wp),                DIMENSION(jpi,jpj,1)         :: tmask_2d 
     136      REAL(wp),                DIMENSION(jpi,jpj)           :: pmld_2d 
     137      REAL(wp),                DIMENSION(jpi,jpj)           :: mld_max_bkg_2d 
    153138      !!--------------------------------------------------------------------------- 
    154139 
    155140      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value 
    156141      IF ( p_maxchlinc > 0.0 ) THEN 
    157          IF ( ld_chltot ) THEN 
    158             DO jj = 1, jpj 
    159                DO ji = 1, jpi 
    160                   pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) ) 
    161                END DO 
    162             END DO 
    163          ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN 
    164             DO jj = 1, jpj 
    165                DO ji = 1, jpi 
    166                   pinc_chltot(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) 
    167                   pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) ) 
    168                   IF ( pinc_chltot(ji,jj) .NE. ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) ) ) THEN 
    169                      zfrac = pinc_chltot(ji,jj) / ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) ) 
    170                      pinc_chldia(ji,jj) = pinc_chldia(ji,jj) * zfrac 
    171                      pinc_chlnon(ji,jj) = pinc_chlnon(ji,jj) * zfrac 
    172                   ENDIF 
    173                END DO 
    174             END DO 
    175          ELSE IF ( ld_chldia ) THEN 
    176             DO jj = 1, jpj 
    177                DO ji = 1, jpi 
    178                   pinc_chldia(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chldia(ji,jj), p_maxchlinc ) ) 
    179                   pinc_chltot(ji,jj) = pinc_chldia(ji,jj) 
    180                END DO 
    181             END DO 
    182          ELSE IF ( ld_chlnon ) THEN 
    183             DO jj = 1, jpj 
    184                DO ji = 1, jpi 
    185                   pinc_chlnon(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chlnon(ji,jj), p_maxchlinc ) ) 
    186                   pinc_chltot(ji,jj) = pinc_chlnon(ji,jj) 
    187                END DO 
    188             END DO 
    189          ENDIF 
    190       ENDIF 
    191  
    192       IF ( ld_phytot .OR. ld_phydia .OR. ld_phynon ) THEN 
    193          CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' ) 
    194       ENDIF 
    195        
    196       IF ( ld_phytobal ) THEN   ! Nitrogen balancing 
    197  
    198          ! Set up model parameters to be passed into Hemmings balancing routine. 
    199          ! For now these are hardwired to the standard HadOCC parameter values 
    200          ! (except C:N ratios) as this is what the scheme was developed for. 
    201          ! Obviously, HadOCC and MEDUSA are rather different models, so this 
    202          ! isn't ideal, but there's not always direct analogues between the two 
    203          ! parameter sets, so it's the easiest way to get something running. 
    204          ! In the longer term, some serious MarMOT-based development is required. 
    205          modparm(1)  = 0.1                               ! grow_sat 
    206          modparm(2)  = 2.0                               ! psmax 
    207          modparm(3)  = 0.845                             ! par 
    208          modparm(4)  = 0.02                              ! alpha 
    209          modparm(5)  = 0.05                              ! resp_rate 
    210          modparm(6)  = 0.05                              ! pmort_rate 
    211          modparm(7)  = 0.01                              ! phyto_min 
    212          modparm(8)  = 0.05                              ! z_mort_1 
    213          modparm(9)  = 1.0                               ! z_mort_2 
    214          modparm(10) = ( xthetapn  + xthetapd  ) / 2.0   ! c2n_p 
    215          modparm(11) = ( xthetazmi + xthetazme ) / 2.0   ! c2n_z 
    216          modparm(12) = xthetad                           ! c2n_d 
    217          modparm(13) = 0.01                              ! graze_threshold 
    218          modparm(14) = 2.0                               ! holling_coef 
    219          modparm(15) = 0.5                               ! graze_sat 
    220          modparm(16) = 2.0                               ! graze_max 
    221  
    222          ! Set up assimilation parameters to be passed into balancing routine 
    223          ! Not sure what assimparm(1) is meant to be, but it doesn't get used 
    224          assimparm(2)  = balnutext 
    225          assimparm(3)  = balnutmin 
    226          assimparm(4)  = r 
    227          assimparm(5)  = beta_g 
    228          assimparm(6)  = beta_l 
    229          assimparm(7)  = beta_m 
    230          assimparm(8)  = a_g 
    231          assimparm(9)  = a_l 
    232          assimparm(10) = a_m 
    233          assimparm(11) = zfracb0 
    234          assimparm(12) = zfracb1 
    235          assimparm(13) = qrfmax 
    236          assimparm(14) = qafmax 
    237          assimparm(15) = zrfmax 
    238          assimparm(16) = zafmax 
    239          assimparm(17) = prfmax 
    240          assimparm(18) = incphymin 
    241          assimparm(19) = integnstep 
    242          assimparm(20) = pthreshold 
    243  
    244          ! Set up external tracer indices array bstate 
    245          i_tracer(1) = 1   ! nutrient 
    246          i_tracer(2) = 2   ! phytoplankton 
    247          i_tracer(3) = 3   ! zooplankton 
    248          i_tracer(4) = 4   ! detritus 
    249          i_tracer(5) = 5   ! DIC 
    250          i_tracer(6) = 6   ! Alkalinity 
    251  
    252          ! Set background state 
    253          bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jpdin) 
    254          bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jpphn) + tracer_bkg(:,:,:,jpphd) 
    255          bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jpzmi) + tracer_bkg(:,:,:,jpzme) 
    256          bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jpdet) 
    257          bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jpdic) 
    258          bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jpalk) 
    259  
    260          ! Calculate carbon to chlorophyll ratio for combined phytoplankton 
    261          ! and nitrogen to biomass equivalent for PZD 
    262          ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA 
    263          cchl_p(:,:) = 0.0 
    264          DO jj = 1, jpj 
    265             DO ji = 1, jpi 
    266                IF ( ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) ) .GT. 0.0 ) THEN 
    267                   cchl_p(ji,jj) = xmassc * ( ( tracer_bkg(ji,jj,1,jpphn) * xthetapn ) +      & 
    268                      &                       ( tracer_bkg(ji,jj,1,jpphd) * xthetapd )   ) /  & 
    269                      &            ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) ) 
    270                ENDIF 
    271             END DO 
    272          END DO 
    273          n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) ) 
    274          n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) ) 
    275          n2be_d = 14.01 + ( xmassc * xthetad ) 
    276  
    277          ! Call nitrogen balancing routine 
    278          CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   & 
    279             &               n2be_p, n2be_z, n2be_d, assimparm,                      & 
    280             &               INT(pincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       & 
    281             &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot(:,:), cchl_p(:,:), & 
    282             &               nbal_active, phyt_avg_bkg(:,:),                         & 
    283             &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      & 
    284             &               subsurf_active, deepneg_active,                         & 
    285             &               deeppos_active, nutprof_active,                         & 
    286             &               bstate, outincs,                                        & 
    287             &               diag_active, diag,                                      & 
    288             &               diag_fulldepth_active, diag_fulldepth ) 
    289           
    290          ! Loop over each grid point partioning the increments 
    291          phyto2d_balinc(:,:,:,:) = 0.0 
    292142         DO jk = 1, jpk 
    293143            DO jj = 1, jpj 
    294144               DO ji = 1, jpi 
    295  
    296                   ! Phytoplankton 
    297                   IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. & 
    298                      & ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) .AND. & 
    299                      & ( pinc_chltot(ji,jj) /= 0.0 ) ) THEN 
    300                      IF ( ld_chltot ) THEN 
    301                         ! Phytoplankton nitrogen split up based on existing ratios 
    302                         zfrac_phn = tracer_bkg(ji,jj,jk,jpphn) / & 
    303                            &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd)) 
    304                         zfrac_phd = tracer_bkg(ji,jj,jk,jpphd) / & 
    305                            &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd)) 
    306                      ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN 
    307                         ! Phytoplankton nitrogen split up based on assimilation increments 
    308                         zfrac_phn = pinc_chlnon(ji,jj) / pinc_chltot(ji,jj) 
    309                         zfrac_phd = pinc_chldia(ji,jj) / pinc_chltot(ji,jj) 
    310                      ENDIF 
    311  
    312                      ! Phytoplankton silicate split up based on existing ratios 
    313                      zrat_pds_phd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpphd) 
    314                       
    315                      ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen 
    316                      ! Not using pinc_chltot directly as it's only 2D 
    317                      ! This method should give same results at surface as splitting pinc_chltot would 
    318                      zrat_chn_phn = tracer_bkg(ji,jj,jk,jpchn) / tracer_bkg(ji,jj,jk,jpphn) 
    319                      zrat_chd_phd = tracer_bkg(ji,jj,jk,jpchd) / tracer_bkg(ji,jj,jk,jpphd) 
    320                       
    321                      phyto2d_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn 
    322                      phyto2d_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd 
    323                      phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_pds_phd 
    324                      phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,jk,jpphn) * zrat_chn_phn 
    325                      phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_chd_phd 
    326                   ENDIF 
    327  
    328                   ! Zooplankton nitrogen split up based on existing ratios 
    329                   IF ( ( tracer_bkg(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpzme) > 0.0 ) ) THEN 
    330                      zfrac_zmi = tracer_bkg(ji,jj,jk,jpzmi) / & 
    331                         &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme)) 
    332                      zfrac_zme = tracer_bkg(ji,jj,jk,jpzme) / & 
    333                         &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme)) 
    334                      phyto2d_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi 
    335                      phyto2d_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme 
    336                   ENDIF 
    337  
    338                   ! Nitrogen nutrient straight from balancing scheme 
    339                   phyto2d_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1)) 
    340  
    341                   ! Nitrogen detritus straight from balancing scheme 
    342                   phyto2d_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4)) 
    343  
    344                   ! DIC straight from balancing scheme 
    345                   phyto2d_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5)) 
    346  
    347                   ! Alkalinity straight from balancing scheme 
    348                   phyto2d_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6)) 
    349  
    350                   ! Remove diatom silicate increment from nutrient silicate to conserve mass 
    351                   IF ( ( tracer_bkg(ji,jj,jk,jpsil) - phyto2d_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN 
    352                      phyto2d_balinc(ji,jj,jk,jpsil) = phyto2d_balinc(ji,jj,jk,jppds) * (-1.0) 
    353                   ENDIF 
    354  
    355                   ! Carbon detritus based on existing ratios 
    356                   IF ( ( tracer_bkg(ji,jj,jk,jpdet) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpdtc) > 0.0 ) ) THEN 
    357                      zrat_dtc_det = tracer_bkg(ji,jj,jk,jpdtc) / tracer_bkg(ji,jj,jk,jpdet) 
    358                      phyto2d_balinc(ji,jj,jk,jpdtc) = phyto2d_balinc(ji,jj,jk,jpdet) * zrat_dtc_det 
    359                   ENDIF 
    360  
    361                   ! Do nothing with iron or oxygen for the time being 
    362                   phyto2d_balinc(ji,jj,jk,jpfer) = 0.0 
    363                   phyto2d_balinc(ji,jj,jk,jpoxy) = 0.0 
    364                    
     145                  pinc_chltot(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj,jk), p_maxchlinc ) ) 
    365146               END DO 
    366147            END DO 
    367148         END DO 
    368        
    369       ELSE   ! No nitrogen balancing 
    370        
    371          ! Initialise individual chlorophyll increments to zero 
    372          phyto2d_balinc(:,:,:,jpchn) = 0.0 
    373          phyto2d_balinc(:,:,:,jpchd) = 0.0 
    374           
    375          ! Split up total surface chlorophyll increments 
     149      ENDIF 
     150 
     151      ! Set up model parameters to be passed into Hemmings balancing routine. 
     152      ! For now these are hardwired to the standard HadOCC parameter values 
     153      ! (except C:N ratios) as this is what the scheme was developed for. 
     154      ! Obviously, HadOCC and MEDUSA are rather different models, so this 
     155      ! isn't ideal, but there's not always direct analogues between the two 
     156      ! parameter sets, so it's the easiest way to get something running. 
     157      ! In the longer term, some serious MarMOT-based development is required. 
     158      modparm(1)  = 0.1                               ! grow_sat 
     159      modparm(2)  = 2.0                               ! psmax 
     160      modparm(3)  = 0.845                             ! par 
     161      modparm(4)  = 0.02                              ! alpha 
     162      modparm(5)  = 0.05                              ! resp_rate 
     163      modparm(6)  = 0.05                              ! pmort_rate 
     164      modparm(7)  = 0.01                              ! phyto_min 
     165      modparm(8)  = 0.05                              ! z_mort_1 
     166      modparm(9)  = 1.0                               ! z_mort_2 
     167      modparm(10) = ( xthetapn  + xthetapd  ) / 2.0   ! c2n_p 
     168      modparm(11) = ( xthetazmi + xthetazme ) / 2.0   ! c2n_z 
     169      modparm(12) = xthetad                           ! c2n_d 
     170      modparm(13) = 0.01                              ! graze_threshold 
     171      modparm(14) = 2.0                               ! holling_coef 
     172      modparm(15) = 0.5                               ! graze_sat 
     173      modparm(16) = 2.0                               ! graze_max 
     174 
     175      ! Set up assimilation parameters to be passed into balancing routine 
     176      ! Not sure what assimparm(1) is meant to be, but it doesn't get used 
     177      assimparm(2)  = balnutext 
     178      assimparm(3)  = balnutmin 
     179      assimparm(4)  = r 
     180      assimparm(5)  = beta_g 
     181      assimparm(6)  = beta_l 
     182      assimparm(7)  = beta_m 
     183      assimparm(8)  = a_g 
     184      assimparm(9)  = a_l 
     185      assimparm(10) = a_m 
     186      assimparm(11) = zfracb0 
     187      assimparm(12) = zfracb1 
     188      assimparm(13) = qrfmax 
     189      assimparm(14) = qafmax 
     190      assimparm(15) = zrfmax 
     191      assimparm(16) = zafmax 
     192      assimparm(17) = prfmax 
     193      assimparm(18) = incphymin 
     194      assimparm(19) = integnstep 
     195      assimparm(20) = pthreshold 
     196 
     197      ! Set up external tracer indices array bstate 
     198      i_tracer(1) = 1   ! nutrient 
     199      i_tracer(2) = 2   ! phytoplankton 
     200      i_tracer(3) = 3   ! zooplankton 
     201      i_tracer(4) = 4   ! detritus 
     202      i_tracer(5) = 5   ! DIC 
     203      i_tracer(6) = 6   ! Alkalinity 
     204 
     205      ! Set background state 
     206      bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jpdin) 
     207      bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jpphn) + tracer_bkg(:,:,:,jpphd) 
     208      bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jpzmi) + tracer_bkg(:,:,:,jpzme) 
     209      bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jpdet) 
     210      bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jpdic) 
     211      bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jpalk) 
     212 
     213      ! Calculate carbon to chlorophyll ratio for combined phytoplankton 
     214      ! and nitrogen to biomass equivalent for PZD 
     215      ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA 
     216      cchl_p(:,:,:) = 0.0 
     217      DO jk = 1, jpk 
    376218         DO jj = 1, jpj 
    377219            DO ji = 1, jpi 
    378                IF ( ( tracer_bkg(ji,jj,1,jpchn) > 0.0 ) .AND. & 
    379                   & ( tracer_bkg(ji,jj,1,jpchd) > 0.0 ) ) THEN 
    380                   IF ( ld_chltot ) THEN 
    381                      ! Chlorophyll split up based on existing ratios 
    382                      zfrac_chn = tracer_bkg(ji,jj,1,jpchn) / & 
    383                         &        ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) ) 
    384                      zfrac_chd = tracer_bkg(ji,jj,1,jpchd) / & 
    385                         &        ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) ) 
    386                      phyto2d_balinc(ji,jj,1,jpchn) = pinc_chltot(ji,jj) * zfrac_chn 
    387                      phyto2d_balinc(ji,jj,1,jpchd) = pinc_chltot(ji,jj) * zfrac_chd 
    388                   ENDIF 
    389                   IF( ld_chldia ) THEN 
    390                      phyto2d_balinc(ji,jj,1,jpchd) = pinc_chldia(ji,jj) 
    391                   ENDIF 
    392                   IF( ld_chlnon ) THEN 
    393                      phyto2d_balinc(ji,jj,1,jpchn) = pinc_chlnon(ji,jj) 
    394                   ENDIF 
    395                    
    396                   ! Maintain stoichiometric ratios of nitrogen and silicate 
    397                   IF ( ld_chltot .OR. ld_chlnon ) THEN 
    398                      zrat_phn_chn = tracer_bkg(ji,jj,1,jpphn) / tracer_bkg(ji,jj,1,jpchn) 
    399                      phyto2d_balinc(ji,jj,1,jpphn) = phyto2d_balinc(ji,jj,1,jpchn) * zrat_phn_chn 
    400                   ENDIF 
    401                   IF ( ld_chltot .OR. ld_chldia ) THEN 
    402                      zrat_phd_chd = tracer_bkg(ji,jj,1,jpphd) / tracer_bkg(ji,jj,1,jpchd) 
    403                      phyto2d_balinc(ji,jj,1,jpphd) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_phd_chd 
    404                      zrat_pds_chd = tracer_bkg(ji,jj,1,jppds) / tracer_bkg(ji,jj,1,jpchd) 
    405                      phyto2d_balinc(ji,jj,1,jppds) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_pds_chd 
    406                   ENDIF 
     220               IF ( ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd ) ) .GT. 0.0 ) THEN 
     221                  cchl_p(ji,jj,jk) = xmassc * ( ( tracer_bkg(ji,jj,jk,jpphn) * xthetapn ) +      & 
     222                     &                       ( tracer_bkg(ji,jj,jk,jpphd) * xthetapd )   ) /  & 
     223                     &            ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd ) ) 
    407224               ENDIF 
    408225            END DO 
    409226         END DO 
     227      END DO 
     228      n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) ) 
     229      n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) ) 
     230      n2be_d = 14.01 + ( xmassc * xthetad ) 
     231 
     232      pmld_2d(:,:)        = 0.5 
     233      mld_max_bkg_2d(:,:) = 0.5 
     234 
     235      DO jk = 1, jpk 
     236 
     237         tmask_2d(:,:,1)       = tmask(:,:,jk) 
     238         pinc_chltot_2d(:,:)   = pinc_chltot(:,:,jk) 
     239         cchl_p_2d(:,:)        = cchl_p(:,:,jk) 
     240         phyt_avg_bkg_2d(:,:)  = phyt_avg_bkg(:,:,jk) 
     241         pgrow_avg_bkg_2d(:,:) = pgrow_avg_bkg(:,:,jk) 
     242         ploss_avg_bkg_2d(:,:) = ploss_avg_bkg(:,:,jk) 
     243         bstate_2d(:,:,1,:)    = bstate(:,:,jk,:) 
     244         outincs_2d(:,:,:,:)   = 0.0 
    410245          
    411          ! Propagate through mixed layer 
     246         ! Call nitrogen balancing routine 
     247         CALL bio_analysis( jpi, jpj, 1, gdepw_n(:,:,2), i_tracer, modparm,   & 
     248            &               n2be_p, n2be_z, n2be_d, assimparm,                      & 
     249            &               INT(pincper), 1, INT(SUM(tmask_2d,3)), tmask_2d(:,:,:),       & 
     250            &               pmld_2d(:,:), mld_max_bkg_2d(:,:), pinc_chltot_2d(:,:), cchl_p_2d(:,:), & 
     251            &               nbal_active, phyt_avg_bkg_2d(:,:),                         & 
     252            &               gl_active, pgrow_avg_bkg_2d(:,:), ploss_avg_bkg_2d(:,:),      & 
     253            &               subsurf_active, deepneg_active,                         & 
     254            &               deeppos_active, nutprof_active,                         & 
     255            &               bstate_2d, outincs_2d,                                        & 
     256            &               diag_active, diag,                                      & 
     257            &               diag_fulldepth_active, diag_fulldepth ) 
     258 
     259         outincs(:,:,jk,:) = outincs_2d(:,:,1,:) 
     260 
     261      END DO 
     262 
     263      ! Loop over each grid point partioning the increments 
     264      phyto3d_balinc(:,:,:,:) = 0.0 
     265      DO jk = 1, jpk 
    412266         DO jj = 1, jpj 
    413267            DO ji = 1, jpi 
    414                ! 
    415                jkmax = jpk-1 
    416                DO jk = jpk-1, 1, -1 
    417                   IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. & 
    418                      & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN 
    419                      pmld(ji,jj) = gdepw_n(ji,jj,jk+1) 
    420                      jkmax = jk 
    421                   ENDIF 
    422                END DO 
    423                ! 
    424                DO jk = 2, jkmax 
    425                   phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,1,jpchn) 
    426                   phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,1,jpchd) 
    427                   phyto2d_balinc(ji,jj,jk,jpphn) = phyto2d_balinc(ji,jj,1,jpphn) 
    428                   phyto2d_balinc(ji,jj,jk,jpphd) = phyto2d_balinc(ji,jj,1,jpphd) 
    429                   phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,1,jppds) 
    430                END DO 
    431                ! 
     268 
     269               ! Phytoplankton 
     270               IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. & 
     271                  & ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) .AND. & 
     272                  & ( pinc_chltot(ji,jj,jk) /= 0.0 ) ) THEN 
     273                  ! Phytoplankton nitrogen split up based on existing ratios 
     274                  zfrac_phn = tracer_bkg(ji,jj,jk,jpphn) / & 
     275                     &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd)) 
     276                  zfrac_phd = tracer_bkg(ji,jj,jk,jpphd) / & 
     277                     &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd)) 
     278 
     279                  ! Phytoplankton silicate split up based on existing ratios 
     280                  zrat_pds_phd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpphd) 
     281 
     282                  ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen 
     283                  ! Not using pinc_chltot directly as it's only 2D 
     284                  ! This method should give same results at surface as splitting pinc_chltot would 
     285                  zrat_chn_phn = tracer_bkg(ji,jj,jk,jpchn) / tracer_bkg(ji,jj,jk,jpphn) 
     286                  zrat_chd_phd = tracer_bkg(ji,jj,jk,jpchd) / tracer_bkg(ji,jj,jk,jpphd) 
     287 
     288                  phyto3d_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn 
     289                  phyto3d_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd 
     290                  phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpphd) * zrat_pds_phd 
     291                  phyto3d_balinc(ji,jj,jk,jpchn) = phyto3d_balinc(ji,jj,jk,jpphn) * zrat_chn_phn 
     292                  phyto3d_balinc(ji,jj,jk,jpchd) = phyto3d_balinc(ji,jj,jk,jpphd) * zrat_chd_phd 
     293               ENDIF 
     294 
     295               ! Zooplankton nitrogen split up based on existing ratios 
     296               IF ( ( tracer_bkg(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpzme) > 0.0 ) ) THEN 
     297                  zfrac_zmi = tracer_bkg(ji,jj,jk,jpzmi) / & 
     298                     &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme)) 
     299                  zfrac_zme = tracer_bkg(ji,jj,jk,jpzme) / & 
     300                     &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme)) 
     301                  phyto3d_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi 
     302                  phyto3d_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme 
     303               ENDIF 
     304 
     305               ! Nitrogen nutrient straight from balancing scheme 
     306               phyto3d_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1)) 
     307 
     308               ! Nitrogen detritus straight from balancing scheme 
     309               phyto3d_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4)) 
     310 
     311               ! DIC straight from balancing scheme 
     312               phyto3d_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5)) 
     313 
     314               ! Alkalinity straight from balancing scheme 
     315               phyto3d_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6)) 
     316 
     317               ! Remove diatom silicate increment from nutrient silicate to conserve mass 
     318               IF ( ( tracer_bkg(ji,jj,jk,jpsil) - phyto3d_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN 
     319                  phyto3d_balinc(ji,jj,jk,jpsil) = phyto3d_balinc(ji,jj,jk,jppds) * (-1.0) 
     320               ENDIF 
     321 
     322               ! Carbon detritus based on existing ratios 
     323               IF ( ( tracer_bkg(ji,jj,jk,jpdet) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpdtc) > 0.0 ) ) THEN 
     324                  zrat_dtc_det = tracer_bkg(ji,jj,jk,jpdtc) / tracer_bkg(ji,jj,jk,jpdet) 
     325                  phyto3d_balinc(ji,jj,jk,jpdtc) = phyto3d_balinc(ji,jj,jk,jpdet) * zrat_dtc_det 
     326               ENDIF 
     327 
     328               ! Do nothing with iron or oxygen for the time being 
     329               phyto3d_balinc(ji,jj,jk,jpfer) = 0.0 
     330               phyto3d_balinc(ji,jj,jk,jpoxy) = 0.0 
     331 
    432332            END DO 
    433333         END DO 
    434  
    435          ! Set other balancing increments to zero 
    436          phyto2d_balinc(:,:,:,jpzmi) = 0.0 
    437          phyto2d_balinc(:,:,:,jpzme) = 0.0 
    438          phyto2d_balinc(:,:,:,jpdin) = 0.0 
    439          phyto2d_balinc(:,:,:,jpsil) = 0.0 
    440          phyto2d_balinc(:,:,:,jpfer) = 0.0 
    441          phyto2d_balinc(:,:,:,jpdet) = 0.0 
    442          phyto2d_balinc(:,:,:,jpdtc) = 0.0 
    443          phyto2d_balinc(:,:,:,jpdic) = 0.0 
    444          phyto2d_balinc(:,:,:,jpalk) = 0.0 
    445          phyto2d_balinc(:,:,:,jpoxy) = 0.0 
    446  
    447       ENDIF 
     334      END DO 
    448335       
    449336      ! If performing extra tidal mixing in the Indonesian Throughflow, 
     
    453340         DO jn = 1, jptra 
    454341            DO jk = 1, jpk 
    455                phyto2d_balinc(:,:,jk,jn) = phyto2d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) ) 
     342               phyto3d_balinc(:,:,jk,jn) = phyto3d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) ) 
    456343            END DO 
    457344         END DO 
    458345      ENDIF 
    459346 
    460    END SUBROUTINE asm_phyto2d_bal_medusa 
     347   END SUBROUTINE asm_phyto3d_bal_medusa 
    461348 
    462349#else 
     
    465352   !!---------------------------------------------------------------------- 
    466353CONTAINS 
    467    SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      & 
    468       &                              pinc_chltot,                    & 
    469       &                              ld_chldia,                      & 
    470       &                              pinc_chldia,                    & 
    471       &                              ld_chlnon,                      & 
    472       &                              pinc_chlnon,                    & 
    473       &                              ld_phytot,                      & 
    474       &                              pinc_phytot,                    & 
    475       &                              ld_phydia,                      & 
    476       &                              pinc_phydia,                    & 
    477       &                              ld_phynon,                      & 
    478       &                              pinc_phynon,                    & 
     354   SUBROUTINE asm_phyto3d_bal_medusa(pinc_chltot,                    & 
    479355      &                              pincper,                        & 
    480       &                              p_maxchlinc, ld_phytobal, pmld, & 
     356      &                              p_maxchlinc,                   & 
    481357      &                              pgrow_avg_bkg, ploss_avg_bkg,   & 
    482       &                              phyt_avg_bkg, mld_max_bkg,      & 
    483       &                              tracer_bkg, phyto2d_balinc ) 
    484       LOGICAL :: ld_chltot 
    485       REAL    :: pinc_chltot(:,:) 
    486       LOGICAL :: ld_chldia 
    487       REAL    :: pinc_chldia(:,:) 
    488       LOGICAL :: ld_chlnon 
    489       REAL    :: pinc_chlnon(:,:) 
    490       LOGICAL :: ld_phytot 
    491       REAL    :: pinc_phytot(:,:) 
    492       LOGICAL :: ld_phydia 
    493       REAL    :: pinc_phydia(:,:) 
    494       LOGICAL :: ld_phynon 
    495       REAL    :: pinc_phynon(:,:) 
     358      &                              phyt_avg_bkg,                   & 
     359      &                              tracer_bkg, phyto3d_balinc ) 
     360      REAL    :: pinc_chltot(:,:,:) 
    496361      REAL    :: pincper 
    497362      REAL    :: p_maxchlinc 
    498       LOGICAL :: ld_phytobal 
    499       REAL    :: pmld(:,:) 
    500       REAL    :: pgrow_avg_bkg(:,:) 
    501       REAL    :: ploss_avg_bkg(:,:) 
    502       REAL    :: phyt_avg_bkg(:,:) 
    503       REAL    :: mld_max_bkg(:,:) 
     363      REAL    :: pgrow_avg_bkg(:,:,:) 
     364      REAL    :: ploss_avg_bkg(:,:,:) 
     365      REAL    :: phyt_avg_bkg(:,:,:) 
    504366      REAL    :: tracer_bkg(:,:,:,:) 
    505367      REAL    :: phyto2d_balinc(:,:,:,:) 
    506       WRITE(*,*) 'asm_phyto2d_bal_medusa: You should not have seen this print! error?' 
    507    END SUBROUTINE asm_phyto2d_bal_medusa 
     368      WRITE(*,*) 'asm_phyto3d_bal_medusa: You should not have seen this print! error?' 
     369   END SUBROUTINE asm_phyto3d_bal_medusa 
    508370#endif 
    509371 
    510372   !!====================================================================== 
    511 END MODULE asmphyto2dbal_medusa 
     373END MODULE asmphyto3dbal_medusa 
Note: See TracChangeset for help on using the changeset viewer.