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 8436 for branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC – NEMO

Ignore:
Timestamp:
2017-08-14T15:22:09+02:00 (7 years ago)
Author:
dford
Message:

Implement initial version of surface chlorophyll assimilation for MEDUSA.

Location:
branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM
Files:
3 edited
1 copied

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90

    r8400 r8436  
    4949#if defined key_lim3 
    5050   USE ice 
     51#endif 
     52#if defined key_hadocc 
     53   USE trc, ONLY: trn, & 
     54      &     pgrow_avg, & 
     55      &     ploss_avg, & 
     56      &     phyt_avg,  & 
     57      &     mld_max,   & 
     58      &     HADOCC_CHL 
     59   USE had_bgc_const, ONLY: cchl_p 
     60   USE par_hadocc 
     61#elif defined key_medusa && defined key_foam_medusa 
     62   USE trc, ONLY: trn 
     63   USE sms_medusa, ONLY: pgrow_avg, & 
     64      &                  ploss_avg, & 
     65      &                  phyt_avg,  & 
     66      &                  mld_max 
     67   USE par_medusa 
    5168#endif 
    5269   IMPLICIT NONE 
     
    121138!            CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               ) 
    122139            CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               ) 
     140#if defined key_hadocc 
     141            CALL iom_rstput( kt, nitbkg_r, inum, 'pgrow_avg'     , pgrow_avg             ) 
     142            CALL iom_rstput( kt, nitbkg_r, inum, 'ploss_avg'     , ploss_avg             ) 
     143            CALL iom_rstput( kt, nitbkg_r, inum, 'phyt_avg'      , phyt_avg              ) 
     144            CALL iom_rstput( kt, nitbkg_r, inum, 'mld_max'       , mld_max               ) 
     145            CALL iom_rstput( kt, nitbkg_r, inum, 'nutrients'     , trn(:,:,:,jp_had_nut) ) 
     146            CALL iom_rstput( kt, nitbkg_r, inum, 'phytoplankton' , trn(:,:,:,jp_had_phy) ) 
     147            CALL iom_rstput( kt, nitbkg_r, inum, 'zooplankton'   , trn(:,:,:,jp_had_zoo) ) 
     148            CALL iom_rstput( kt, nitbkg_r, inum, 'detritus'      , trn(:,:,:,jp_had_pdn) ) 
     149            CALL iom_rstput( kt, nitbkg_r, inum, 'dic'           , trn(:,:,:,jp_had_dic) ) 
     150            CALL iom_rstput( kt, nitbkg_r, inum, 'alkalinity'    , trn(:,:,:,jp_had_alk) ) 
     151            CALL iom_rstput( kt, nitbkg_r, inum, 'chlorophyll'   , HADOCC_CHL(:,:,1)     ) 
     152            CALL iom_rstput( kt, nitbkg_r, inum, 'c_to_chl'      , cchl_p(:,:,1)         ) 
     153#elif defined key_medusa && defined key_foam_medusa 
     154            CALL iom_rstput( kt, nitbkg_r, inum, 'pgrow_avg'     , pgrow_avg                           ) 
     155            CALL iom_rstput( kt, nitbkg_r, inum, 'ploss_avg'     , ploss_avg                           ) 
     156            CALL iom_rstput( kt, nitbkg_r, inum, 'phyt_avg'      , phyt_avg                            ) 
     157            CALL iom_rstput( kt, nitbkg_r, inum, 'mld_max'       , mld_max                             ) 
     158            CALL iom_rstput( kt, nitbkg_r, inum, 'nutrients'     , trn(:,:,:,jpdin)                    ) 
     159            CALL iom_rstput( kt, nitbkg_r, inum, 'phytoplankton' , trn(:,:,:,jpphn) + trn(:,:,:,jpphd) ) 
     160            CALL iom_rstput( kt, nitbkg_r, inum, 'zooplankton'   , trn(:,:,:,jpzmi) + trn(:,:,:,jpzme) ) 
     161            CALL iom_rstput( kt, nitbkg_r, inum, 'detritus'      , trn(:,:,:,jpdet)                    ) 
     162            CALL iom_rstput( kt, nitbkg_r, inum, 'dic'           , trn(:,:,:,jpdic)                    ) 
     163            CALL iom_rstput( kt, nitbkg_r, inum, 'alkalinity'    , trn(:,:,:,jpalk)                    ) 
     164            CALL iom_rstput( kt, nitbkg_r, inum, 'chlorophyll'   , trn(:,:,1,jpchn) + trn(:,:,1,jpchd) ) 
     165#endif 
    123166            ! 
    124167            CALL iom_close( inum ) 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r8428 r8436  
    126126   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: logchl_balinc  !: Increment to BGC variables from logchl assim 
    127127#endif 
     128#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa) 
     129   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pgrow_avg_bkg  !: Background phyto growth for logchl balancing 
     130   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ploss_avg_bkg  !: Background phyto loss for logchl balancing 
     131   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: phyt_avg_bkg   !: Background phyto for logchl balancing 
     132   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mld_max_bkg    !: Background max MLD for logchl balancing 
     133#endif 
    128134   REAL(wp) :: rn_maxchlinc = -999.0  !: maximum absolute non-log chlorophyll increment from logchl assimilation 
    129135                                      !: <= 0 implies no maximum applied (switch turned off) 
     
    620626 
    621627      ENDIF 
     628 
     629#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa) 
     630      IF ( ln_logchlinc ) THEN 
     631 
     632         ALLOCATE( pgrow_avg_bkg(jpi,jpj) ) 
     633         ALLOCATE( ploss_avg_bkg(jpi,jpj) ) 
     634         ALLOCATE( phyt_avg_bkg(jpi,jpj)  ) 
     635         ALLOCATE( mld_max_bkg(jpi,jpj)   ) 
     636         pgrow_avg_bkg(:,:) = 0.0 
     637         ploss_avg_bkg(:,:) = 0.0 
     638         phyt_avg_bkg(:,:)  = 0.0 
     639         mld_max_bkg(:,:)   = 0.0 
     640          
     641         IF ( ln_logchlbal ) THEN 
     642          
     643            !-------------------------------------------------------------------- 
     644            ! Read background variables for logchl balancing 
     645            !-------------------------------------------------------------------- 
     646 
     647            CALL iom_open( c_asmbkg, inum ) 
     648 
     649            CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg ) 
     650            CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg ) 
     651            CALL iom_get( inum, jpdom_autoglo, 'phyt_avg',  phyt_avg_bkg  ) 
     652            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   ) 
     653            pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1) 
     654            ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1) 
     655            phyt_avg_bkg(:,:)  = phyt_avg_bkg(:,:)  * tmask(:,:,1) 
     656            mld_max_bkg(:,:)   = mld_max_bkg(:,:)   * tmask(:,:,1) 
     657 
     658            CALL iom_close( inum ) 
     659          
     660         ENDIF 
     661       
     662      ENDIF 
     663#endif 
    622664      ! 
    623665   END SUBROUTINE asm_inc_init 
     
    12661308 
    12671309#if defined key_medusa && defined key_foam_medusa 
    1268          !CALL asm_logchl_bal_medusa() 
    1269          CALL ctl_stop( 'Attempting to assimilate logchl into MEDUSA, ', & 
    1270             &           'but not fully implemented yet' ) 
     1310         CALL asm_logchl_bal_medusa( logchl_bkginc, zincper, mld_choice_bgc, & 
     1311            &                        rn_maxchlinc, ln_logchlbal,             & 
     1312            &                        pgrow_avg_bkg, ploss_avg_bkg,           & 
     1313            &                        phyt_avg_bkg, mld_max_bkg,              & 
     1314            &                        logchl_balinc ) 
    12711315#elif defined key_hadocc 
    12721316         CALL asm_logchl_bal_hadocc( logchl_bkginc, zincper, mld_choice_bgc, & 
    1273             &                        rn_maxchlinc, ln_logchlbal, logchl_balinc ) 
     1317            &                        rn_maxchlinc, ln_logchlbal,             & 
     1318            &                        pgrow_avg_bkg, ploss_avg_bkg,           & 
     1319            &                        phyt_avg_bkg, mld_max_bkg,              & 
     1320            &                        logchl_balinc ) 
    12741321#else 
    12751322         CALL ctl_stop( 'Attempting to assimilate logchl, ', & 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_hadocc.F90

    r8428 r8436  
    2525   USE iom                                 ! i/o 
    2626   USE trc,           ONLY: trn, trb,    & ! HadOCC variables 
    27       &                     HADOCC_CHL,  & 
    28       &                     pgrow_avg,   & 
    29       &                     ploss_avg,   & 
    30       &                     phyt_avg,    & 
    31       &                     mld_max 
     27      &                     HADOCC_CHL 
    3228   USE par_hadocc                          ! HadOCC parameters 
    3329   USE had_bgc_stnd,  ONLY: kmt            ! HadOCC parameters 
     
    7571 
    7672   SUBROUTINE asm_logchl_bal_hadocc( logchl_bkginc, aincper, mld_choice_bgc, & 
    77       &                              k_maxchlinc, ld_logchlbal, logchl_balinc ) 
     73      &                              k_maxchlinc, ld_logchlbal,              & 
     74      &                              pgrow_avg_bkg, ploss_avg_bkg,           & 
     75      &                              phyt_avg_bkg, mld_max_bkg,              & 
     76      &                              logchl_balinc ) 
    7877      !!--------------------------------------------------------------------------- 
    7978      !!                    ***  ROUTINE asm_logchl_bal_hadocc  *** 
     
    9594      REAL(wp), INTENT(in   )                               :: k_maxchlinc    ! Max chl increment 
    9695      LOGICAL,  INTENT(in   )                               :: ld_logchlbal   ! Balancing y/n 
     96      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth 
     97      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss 
     98      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto 
     99      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD 
    97100      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc  ! Balancing increments 
    98101      !! 
    99       INTEGER                                               :: ji, jj, jk    ! Loop counters 
     102      INTEGER                                               :: ji, jj, jk, jn ! Loop counters 
    100103      INTEGER                                               :: jkmax          ! Loop index 
    101104      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices 
     
    215218 
    216219         ! Call nitrogen balancing routine 
    217          CALL bio_analysis( jpi, jpj, jpk, ZDZ(:,:,:), i_tracer, modparm,          & 
    218             &               n2be_p, n2be_z, n2be_d, assimparm,                     & 
    219             &               INT(aincper), 1, kmt(:,:), tmask(:,:,:),               & 
    220             &               zmld(:,:), mld_max(:,:), chl_inc(:,:), cchl_p(:,:,1), & 
    221             &               nbal_active, phyt_avg(:,:),                            & 
    222             &               gl_active, pgrow_avg(:,:), ploss_avg(:,:),             & 
    223             &               subsurf_active, deepneg_active,                        & 
    224             &               deeppos_active, nutprof_active,                        & 
    225             &               bstate, outincs,                                       & 
    226             &               diag_active, diag,                                     & 
     220         CALL bio_analysis( jpi, jpj, jpk, ZDZ(:,:,:), i_tracer, modparm,             & 
     221            &               n2be_p, n2be_z, n2be_d, assimparm,                        & 
     222            &               INT(aincper), 1, kmt(:,:), tmask(:,:,:),                  & 
     223            &               zmld(:,:), mld_max_bkg(:,:), chl_inc(:,:), cchl_p(:,:,1), & 
     224            &               nbal_active, phyt_avg_bkg(:,:),                           & 
     225            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),        & 
     226            &               subsurf_active, deepneg_active,                           & 
     227            &               deeppos_active, nutprof_active,                           & 
     228            &               bstate, outincs,                                          & 
     229            &               diag_active, diag,                                        & 
    227230            &               diag_fulldepth_active, diag_fulldepth ) 
    228231 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_medusa.F90

    r8428 r8436  
    1 MODULE asmlogchlbal_hadocc 
     1MODULE asmlogchlbal_medusa 
    22   !!====================================================================== 
    3    !!                       ***  MODULE asmlogchlbal_hadocc  *** 
    4    !! Calculate increments to HadOCC based on surface logchl increments 
     3   !!                       ***  MODULE asmlogchlbal_medusa  *** 
     4   !! Calculate increments to MEDUSA based on surface logchl increments 
    55   !! 
    66   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al. 
     
    1010   !! 
    1111   !!====================================================================== 
    12    !! History :  3.6  ! 2017-08  (D. Ford)     Adapted from bioanal.F90 
    13    !!---------------------------------------------------------------------- 
    14 #if defined key_asminc && defined key_hadocc 
     12   !! History :  3.6  ! 2017-08 (D. Ford)  Adapted from asmlogchlbal_hadocc 
     13   !!---------------------------------------------------------------------- 
     14#if defined key_asminc && defined key_medusa && defined key_foam_medusa 
    1515   !!---------------------------------------------------------------------- 
    1616   !! 'key_asminc'          : assimilation increment interface 
    17    !! 'key_hadocc'          : HadOCC model 
    18    !!---------------------------------------------------------------------- 
    19    !! asm_logchl_bal_hadocc : routine to calculate increments to HadOCC 
     17   !! 'key_medusa'          : MEDUSA model 
     18   !! 'key_foam_medusa'     : MEDUSA extras for FOAM OBS and ASM 
     19   !!---------------------------------------------------------------------- 
     20   !! asm_logchl_bal_medusa : routine to calculate increments to MEDUSA 
    2021   !!---------------------------------------------------------------------- 
    2122   USE par_kind,      ONLY: wp             ! kind parameters 
     
    2425   USE zdfmxl                              ! mixed layer depth 
    2526   USE iom                                 ! i/o 
    26    USE trc,           ONLY: trn, trb,    & ! HadOCC variables 
    27       &                     HADOCC_CHL,  & 
    28       &                     pgrow_avg,   & 
    29       &                     ploss_avg,   & 
    30       &                     phyt_avg,    & 
    31       &                     mld_max 
    32    USE par_hadocc                          ! HadOCC parameters 
    33    USE had_bgc_stnd,  ONLY: kmt            ! HadOCC parameters 
    34    USE had_bgc_const                       ! HadOCC parameters 
     27   USE trc,           ONLY: trn, trb       ! MEDUSA variables 
     28   USE sms_medusa                          ! MEDUSA parameters 
     29   USE par_medusa                          ! MEDUSA parameters 
    3530   USE par_trc,       ONLY: jptra          ! Tracer parameters 
    3631   USE bioanalysis                         ! Nitrogen balancing 
     
    3934   PRIVATE                    
    4035 
    41    PUBLIC asm_logchl_bal_hadocc 
     36   PUBLIC asm_logchl_bal_medusa 
    4237 
    4338   ! Default values for biological assimilation parameters 
     
    7469CONTAINS 
    7570 
    76    SUBROUTINE asm_logchl_bal_hadocc( logchl_bkginc, aincper, mld_choice_bgc, & 
    77       &                              k_maxchlinc, ld_logchlbal, logchl_balinc ) 
     71   SUBROUTINE asm_logchl_bal_medusa( logchl_bkginc, aincper, mld_choice_bgc, & 
     72      &                              k_maxchlinc, ld_logchlbal,              & 
     73      &                              pgrow_avg_bkg, ploss_avg_bkg,           & 
     74      &                              phyt_avg_bkg, mld_max_bkg,              & 
     75      &                              logchl_balinc ) 
    7876      !!--------------------------------------------------------------------------- 
    79       !!                    ***  ROUTINE asm_logchl_bal_hadocc  *** 
    80       !! 
    81       !! ** Purpose :   calculate increments to HadOCC from logchl increments 
     77      !!                    ***  ROUTINE asm_logchl_bal_medusa  *** 
     78      !! 
     79      !! ** Purpose :   calculate increments to MEDUSA from logchl increments 
    8280      !! 
    8381      !! ** Method  :   convert logchl increments to chl increments 
     82      !!                average up MEDUSA to look like HadOCC 
    8483      !!                call nitrogen balancing scheme 
     84      !!                separate back out to MEDUSA 
    8585      !! 
    8686      !! ** Action  :   populate logchl_balinc 
     
    9595      REAL(wp), INTENT(in   )                               :: k_maxchlinc    ! Max chl increment 
    9696      LOGICAL,  INTENT(in   )                               :: ld_logchlbal   ! Balancing y/n 
     97      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth 
     98      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss 
     99      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto 
     100      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD 
    97101      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc  ! Balancing increments 
    98102      !! 
    99       INTEGER                                               :: ji, jj, jk    ! Loop counters 
     103      INTEGER                                               :: ji, jj, jk, jn ! Loop counters 
    100104      INTEGER                                               :: jkmax          ! Loop index 
    101105      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices 
     106      REAL(wp)                                              :: n2be_p         ! N:biomass for total phy 
     107      REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo 
     108      REAL(wp)                                              :: n2be_d         ! N:biomass for detritus 
     109      REAL(wp)                                              :: zfrac_chn      ! Fraction of jpchn 
     110      REAL(wp)                                              :: zfrac_chd      ! Fraction of jpchd 
     111      REAL(wp)                                              :: zfrac_phn      ! Fraction of jpphn 
     112      REAL(wp)                                              :: zfrac_phd      ! Fraction of jpphd 
     113      REAL(wp)                                              :: zfrac_zmi      ! Fraction of jpzmi 
     114      REAL(wp)                                              :: zfrac_zme      ! Fraction of jpzme 
     115      REAL(wp)                                              :: zrat_pds_phd   ! Ratio of jppds:jpphd 
     116      REAL(wp)                                              :: zrat_chd_phd   ! Ratio of jpchd:jpphd 
     117      REAL(wp)                                              :: zrat_chn_phn   ! Ratio of jpchn:jpphn 
     118      REAL(wp)                                              :: zrat_dtc_det   ! Ratio of jpdtc:jpdet 
    102119      REAL(wp),                DIMENSION(jpi,jpj)           :: chl_inc        ! Chlorophyll increments 
     120      REAL(wp),                DIMENSION(jpi,jpj)           :: medusa_chl     ! MEDUSA total chlorophyll 
     121      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy 
    103122      REAL(wp),                DIMENSION(jpi,jpj)           :: zmld           ! Mixed layer depth 
    104123      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters 
     
    117136      ! 3) Subtract background from analysis to get chl incs 
    118137      ! If k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value 
     138      medusa_chl(:,:) = trb(:,:,1,jpchn) + trb(:,:,1,jpchd) 
    119139      DO jj = 1, jpj 
    120140         DO ji = 1, jpi 
    121             IF ( HADOCC_CHL(ji,jj,1) > 0.0 ) THEN 
    122                chl_inc(ji,jj) = 10**( LOG10( HADOCC_CHL(ji,jj,1) ) + logchl_bkginc(ji,jj) ) - HADOCC_CHL(ji,jj,1) 
     141            IF ( medusa_chl(ji,jj) > 0.0 ) THEN 
     142               chl_inc(ji,jj) = 10**( LOG10( medusa_chl(ji,jj) ) + logchl_bkginc(ji,jj) ) - medusa_chl(ji,jj) 
    123143               IF ( k_maxchlinc > 0.0 ) THEN 
    124144                  chl_inc(ji,jj) = MAX( -1.0 * k_maxchlinc, MIN( chl_inc(ji,jj), k_maxchlinc ) ) 
     
    158178      IF ( ld_logchlbal ) THEN   ! Nitrogen balancing 
    159179 
    160          ! Set up model parameters to be passed into Hemmings balancing routine 
    161          modparm(1)  = grow_sat 
    162          modparm(2)  = psmax 
    163          modparm(3)  = par 
    164          modparm(4)  = alpha 
    165          modparm(5)  = resp_rate 
    166          modparm(6)  = pmort_rate 
    167          modparm(7)  = phyto_min 
    168          modparm(8)  = z_mort_1 
    169          modparm(9)  = z_mort_2 
    170          modparm(10) = c2n_p 
    171          modparm(11) = c2n_z 
    172          modparm(12) = c2n_d 
    173          modparm(13) = graze_threshold 
    174          modparm(14) = holling_coef 
    175          modparm(15) = graze_sat 
    176          modparm(16) = graze_max 
     180         ! Set up model parameters to be passed into Hemmings balancing routine. 
     181         ! For now these are hardwired to the standard HadOCC parameter values 
     182         ! (except C:N ratios) as this is what the scheme was developed for. 
     183         ! Obviously, HadOCC and MEDUSA are rather different models, so this 
     184         ! isn't ideal, but there's not always direct analogues between the two 
     185         ! parameter sets, so it's the easiest way to get something running. 
     186         ! In the longer term, some serious MarMOT-based development is required. 
     187         modparm(1)  = 0.1                               ! grow_sat 
     188         modparm(2)  = 2.0                               ! psmax 
     189         modparm(3)  = 0.845                             ! par 
     190         modparm(4)  = 0.02                              ! alpha 
     191         modparm(5)  = 0.05                              ! resp_rate 
     192         modparm(6)  = 0.05                              ! pmort_rate 
     193         modparm(7)  = 0.01                              ! phyto_min 
     194         modparm(8)  = 0.05                              ! z_mort_1 
     195         modparm(9)  = 1.0                               ! z_mort_2 
     196         modparm(10) = ( xthetapn  + xthetapd  ) / 2.0   ! c2n_p 
     197         modparm(11) = ( xthetazmi + xthetazme ) / 2.0   ! c2n_z 
     198         modparm(12) = xthetad                           ! c2n_d 
     199         modparm(13) = 0.01                              ! graze_threshold 
     200         modparm(14) = 2.0                               ! holling_coef 
     201         modparm(15) = 0.5                               ! graze_sat 
     202         modparm(16) = 2.0                               ! graze_max 
    177203 
    178204         ! Set up assimilation parameters to be passed into balancing routine 
     
    207233 
    208234         ! Set background state 
    209          bstate(:,:,:,i_tracer(1)) = trb(:,:,:,jp_had_nut) 
    210          bstate(:,:,:,i_tracer(2)) = trb(:,:,:,jp_had_phy) 
    211          bstate(:,:,:,i_tracer(3)) = trb(:,:,:,jp_had_zoo) 
    212          bstate(:,:,:,i_tracer(4)) = trb(:,:,:,jp_had_pdn) 
    213          bstate(:,:,:,i_tracer(5)) = trb(:,:,:,jp_had_dic) 
    214          bstate(:,:,:,i_tracer(6)) = trb(:,:,:,jp_had_alk) 
     235         bstate(:,:,:,i_tracer(1)) = trb(:,:,:,jpdin) 
     236         bstate(:,:,:,i_tracer(2)) = trb(:,:,:,jpphn) + trb(:,:,:,jpphd) 
     237         bstate(:,:,:,i_tracer(3)) = trb(:,:,:,jpzmi) + trb(:,:,:,jpzme) 
     238         bstate(:,:,:,i_tracer(4)) = trb(:,:,:,jpdet) 
     239         bstate(:,:,:,i_tracer(5)) = trb(:,:,:,jpdic) 
     240         bstate(:,:,:,i_tracer(6)) = trb(:,:,:,jpalk) 
     241 
     242         ! Calculate carbon to chlorophyll ratio for combined phytoplankton 
     243         ! and nitrogen to biomass equivalent for PZD 
     244         ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA 
     245         !cchl_p(:,:) = ( trb(:,:,1,jpchn) + trb(:,:,1,jpchd ) ) / & 
     246         !   &          ( ( trb(:,:,1,jpphn) * xthetapn ) + ( trb(:,:,1,jpphd) * xthetapd ) ) 
     247         cchl_p(:,:) = 0.0 
     248         DO jj = 1, jpj 
     249            DO ji = 1, jpi 
     250               IF ( ( trb(ji,jj,1,jpchn) + trb(ji,jj,1,jpchd ) ) .GT. 0.0 ) THEN 
     251                  cchl_p(ji,jj) = ( ( trb(ji,jj,1,jpphn) * xthetapn ) + ( trb(ji,jj,1,jpphd) * xthetapd ) ) / & 
     252                     &            ( trb(ji,jj,1,jpchn) + trb(ji,jj,1,jpchd ) ) 
     253               ENDIF 
     254            END DO 
     255         END DO 
     256         n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) ) 
     257         n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) ) 
     258         n2be_d = 14.01 + ( xmassc * xthetad ) 
     259          
     260         WRITE(numout,*) 'DAF: nproc, min/max cchl_p, min/max chl_inc = ', nproc, MINVAL(cchl_p), MAXVAL(cchl_p), MINVAL(chl_inc), MAXVAL(chl_inc) 
    215261 
    216262         ! Call nitrogen balancing routine 
    217          CALL bio_analysis( jpi, jpj, jpk, ZDZ(:,:,:), i_tracer, modparm,          & 
    218             &               n2be_p, n2be_z, n2be_d, assimparm,                     & 
    219             &               INT(aincper), 1, kmt(:,:), tmask(:,:,:),               & 
    220             &               zmld(:,:), mld_max(:,:), chl_inc(:,:), cchl_p(:,:,1), & 
    221             &               nbal_active, phyt_avg(:,:),                            & 
    222             &               gl_active, pgrow_avg(:,:), ploss_avg(:,:),             & 
    223             &               subsurf_active, deepneg_active,                        & 
    224             &               deeppos_active, nutprof_active,                        & 
    225             &               bstate, outincs,                                       & 
    226             &               diag_active, diag,                                     & 
     263         CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   & 
     264            &               n2be_p, n2be_z, n2be_d, assimparm,                      & 
     265            &               INT(aincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       & 
     266            &               zmld(:,:), mld_max_bkg(:,:), chl_inc(:,:), cchl_p(:,:), & 
     267            &               nbal_active, phyt_avg_bkg(:,:),                         & 
     268            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      & 
     269            &               subsurf_active, deepneg_active,                         & 
     270            &               deeppos_active, nutprof_active,                         & 
     271            &               bstate, outincs,                                        & 
     272            &               diag_active, diag,                                      & 
    227273            &               diag_fulldepth_active, diag_fulldepth ) 
    228  
    229          ! Save balancing increments 
    230          logchl_balinc(:,:,:,jp_had_nut) = outincs(:,:,:,i_tracer(1)) 
    231          logchl_balinc(:,:,:,jp_had_phy) = outincs(:,:,:,i_tracer(2)) 
    232          logchl_balinc(:,:,:,jp_had_zoo) = outincs(:,:,:,i_tracer(3)) 
    233          logchl_balinc(:,:,:,jp_had_pdn) = outincs(:,:,:,i_tracer(4)) 
    234          logchl_balinc(:,:,:,jp_had_dic) = outincs(:,:,:,i_tracer(5)) 
    235          logchl_balinc(:,:,:,jp_had_alk) = outincs(:,:,:,i_tracer(6)) 
     274          
     275         WRITE(numout,*) 'DAF: nproc, min/max outincs(phy)1,20 = ', nproc, MINVAL(outincs(:,:,1,i_tracer(2))), MAXVAL(outincs(:,:,1,i_tracer(2))), MINVAL(outincs(:,:,20,i_tracer(2))), MAXVAL(outincs(:,:,20,i_tracer(2))) 
     276 
     277         ! Loop over each grid point partioning the increments 
     278         logchl_balinc(:,:,:,:) = 0.0 
     279         DO jk = 1, jpk 
     280            DO jj = 1, jpj 
     281               DO ji = 1, jpi 
     282 
     283                  IF ( ( trb(ji,jj,jk,jpphn) > 0.0 ) .AND. ( trb(ji,jj,jk,jpphd) > 0.0 ) ) THEN 
     284                     ! Phytoplankton nitrogen and silicate split up based on existing ratios 
     285                     zfrac_phn = trb(ji,jj,jk,jpphn) / (trb(ji,jj,jk,jpphn) + trb(ji,jj,jk,jpphd)) 
     286                     zfrac_phd = 1.0 - zfrac_phn 
     287                     zrat_pds_phd = trb(ji,jj,jk,jppds) / trb(ji,jj,jk,jpphd) 
     288                     logchl_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn 
     289                     logchl_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd 
     290                     logchl_balinc(ji,jj,jk,jppds) = logchl_balinc(ji,jj,jk,jpphd) * zrat_pds_phd 
     291 
     292                     ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen 
     293                     ! Not using chl_inc directly as it's only 2D 
     294                     ! This method should give same results at surface as splitting chl_inc would 
     295                     zrat_chn_phn = trb(ji,jj,jk,jpchn) / trb(ji,jj,jk,jpphn) 
     296                     zrat_chd_phd = trb(ji,jj,jk,jpchd) / trb(ji,jj,jk,jpphd) 
     297                     logchl_balinc(ji,jj,jk,jpchn) = logchl_balinc(ji,jj,jk,jpphn) * zrat_chn_phn 
     298                     logchl_balinc(ji,jj,jk,jpchd) = logchl_balinc(ji,jj,jk,jpphd) * zrat_chd_phd 
     299                  ENDIF 
     300 
     301                  IF ( ( trb(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( trb(ji,jj,jk,jpzme) > 0.0 ) ) THEN 
     302                     ! Zooplankton nitrogen split up based on existing ratios 
     303                     zfrac_zmi = trb(ji,jj,jk,jpzmi) / (trb(ji,jj,jk,jpzmi) + trb(ji,jj,jk,jpzme)) 
     304                     zfrac_zme = 1.0 - zfrac_zmi 
     305                     logchl_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi 
     306                     logchl_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme 
     307                  ENDIF 
     308 
     309                  ! Nitrogen nutrient straight from balancing scheme 
     310                  logchl_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1)) 
     311 
     312                  ! Nitrogen detritus straight from balancing scheme 
     313                  logchl_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4)) 
     314 
     315                  ! DIC straight from balancing scheme 
     316                  logchl_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5)) 
     317 
     318                  ! Alkalinity straight from balancing scheme 
     319                  logchl_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6)) 
     320 
     321                  ! Remove diatom silicate increment from nutrient silicate to conserve mass 
     322                  IF ( ( trb(ji,jj,jk,jpsil) - logchl_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN 
     323                     logchl_balinc(ji,jj,jk,jpsil) = logchl_balinc(ji,jj,jk,jppds) * (-1.0) 
     324                  ENDIF 
     325 
     326                  IF ( ( trb(ji,jj,jk,jpdet) > 0.0 ) .AND. ( trb(ji,jj,jk,jpdtc) > 0.0 ) ) THEN 
     327                     ! Carbon detritus based on existing ratios 
     328                     zrat_dtc_det = trb(ji,jj,jk,jpdtc) / trb(ji,jj,jk,jpdet) 
     329                     logchl_balinc(ji,jj,jk,jpdtc) = logchl_balinc(ji,jj,jk,jpdet) * zrat_dtc_det 
     330                  ENDIF 
     331 
     332                  ! Do nothing with iron or oxygen for the time being 
     333                  logchl_balinc(ji,jj,jk,jpfer) = 0.0 
     334                  logchl_balinc(ji,jj,jk,jpoxy) = 0.0 
     335                   
     336               END DO 
     337            END DO 
     338         END DO 
    236339       
    237340      ELSE   ! No nitrogen balancing 
    238341       
    239          ! Initialise phytoplankton increment to zero 
    240          logchl_balinc(:,:,:,jp_had_phy) = 0.0 
     342         ! Initialise individual chlorophyll increments to zero 
     343         logchl_balinc(:,:,:,jpchn) = 0.0 
     344         logchl_balinc(:,:,:,jpchd) = 0.0 
    241345          
    242          ! Convert surface chlorophyll increment to phytoplankton nitrogen 
    243          logchl_balinc(:,:,1,jp_had_phy) = ( cchl_p(:,:,1) / (mw_carbon * c2n_p) ) * chl_inc(:,:) 
     346         ! Split up total surface chlorophyll increments 
     347         DO jj = 1, jpj 
     348            DO ji = 1, jpi 
     349               IF ( medusa_chl(ji,jj) > 0.0 ) THEN 
     350                  zfrac_chn = trb(ji,jj,1,jpchn) / medusa_chl(ji,jj) 
     351                  zfrac_chd = 1.0 - zfrac_chn 
     352                  logchl_balinc(ji,jj,1,jpchn) = chl_inc(ji,jj) * zfrac_chn 
     353                  logchl_balinc(ji,jj,1,jpchd) = chl_inc(ji,jj) * zfrac_chd 
     354               ENDIF 
     355            END DO 
     356         END DO 
    244357          
    245358         ! Propagate through mixed layer 
     
    257370               ! 
    258371               DO jk = 2, jkmax 
    259                   logchl_balinc(ji,jj,jk,jp_had_phy) = logchl_balinc(ji,jj,1,jp_had_phy) 
     372                  logchl_balinc(ji,jj,jk,jpchn) = logchl_balinc(ji,jj,1,jpchn) 
     373                  logchl_balinc(ji,jj,jk,jpchd) = logchl_balinc(ji,jj,1,jpchd) 
    260374               END DO 
    261375               ! 
     
    264378 
    265379         ! Set other balancing increments to zero 
    266          logchl_balinc(:,:,:,jp_had_nut) = 0.0 
    267          logchl_balinc(:,:,:,jp_had_zoo) = 0.0 
    268          logchl_balinc(:,:,:,jp_had_pdn) = 0.0 
    269          logchl_balinc(:,:,:,jp_had_dic) = 0.0 
    270          logchl_balinc(:,:,:,jp_had_alk) = 0.0 
    271        
     380         logchl_balinc(:,:,:,jpphn) = 0.0 
     381         logchl_balinc(:,:,:,jpphd) = 0.0 
     382         logchl_balinc(:,:,:,jppds) = 0.0 
     383         logchl_balinc(:,:,:,jpzmi) = 0.0 
     384         logchl_balinc(:,:,:,jpzme) = 0.0 
     385         logchl_balinc(:,:,:,jpdin) = 0.0 
     386         logchl_balinc(:,:,:,jpsil) = 0.0 
     387         logchl_balinc(:,:,:,jpfer) = 0.0 
     388         logchl_balinc(:,:,:,jpdet) = 0.0 
     389         logchl_balinc(:,:,:,jpdtc) = 0.0 
     390         logchl_balinc(:,:,:,jpdic) = 0.0 
     391         logchl_balinc(:,:,:,jpalk) = 0.0 
     392         logchl_balinc(:,:,:,jpoxy) = 0.0 
     393 
    272394      ENDIF 
    273395 
    274    END SUBROUTINE asm_logchl_bal_hadocc 
     396   END SUBROUTINE asm_logchl_bal_medusa 
    275397 
    276398#else 
     
    279401   !!---------------------------------------------------------------------- 
    280402CONTAINS 
    281    SUBROUTINE asm_logchl_bal_hadocc( logchl_bkginc, aincper, mld_choice_bgc, & 
     403   SUBROUTINE asm_logchl_bal_medusa( logchl_bkginc, aincper, mld_choice_bgc, & 
    282404      &                              k_maxchlinc, logchl_balinc ) 
    283405      REAL    :: logchl_bkginc(:,:) 
     
    286408      REAL    :: k_maxchlinc 
    287409      REAL(   :: logchl_balinc(:,:,:,:) 
    288       WRITE(*,*) 'asm_logchl_bal_hadocc: You should not have seen this print! error?' 
    289    END SUBROUTINE asm_logchl_bal_hadocc 
     410      WRITE(*,*) 'asm_logchl_bal_medusa: You should not have seen this print! error?' 
     411   END SUBROUTINE asm_logchl_bal_medusa 
    290412#endif 
    291413 
    292414   !!====================================================================== 
    293 END MODULE asmlogchlbal_hadocc 
     415END MODULE asmlogchlbal_medusa 
Note: See TracChangeset for help on using the changeset viewer.