LCOV - code coverage report
Current view: top level - configuration/driver - icedrv_init_column.F90 (source / functions) Hit Total Coverage
Test: 200904-015006:19b4ed4c2c:3:base,travis,quick Lines: 786 906 86.75 %
Date: 2020-09-03 20:12:34 Functions: 6 6 100.00 %

          Line data    Source code
       1             : !=========================================================================
       2             : !
       3             : ! Initialization routines for the column package.
       4             : !
       5             : ! author: Elizabeth C. Hunke, LANL
       6             : !
       7             :       module icedrv_init_column
       8             : 
       9             :       use icedrv_kinds
      10             :       use icedrv_domain_size, only: ncat, nx, ncat
      11             :       use icedrv_domain_size, only: max_ntrcr, nblyr, nilyr, nslyr
      12             :       use icedrv_domain_size, only: n_algae, n_zaero, n_doc, n_dic, n_don
      13             :       use icedrv_domain_size, only: n_fed, n_fep, max_nsw, n_bgc, n_aero
      14             :       use icedrv_constants, only: c1, c2, p5, c0, p1
      15             :       use icedrv_constants, only: nu_diag, nu_nml
      16             :       use icepack_intfc, only: icepack_max_don, icepack_max_doc, icepack_max_dic
      17             :       use icepack_intfc, only: icepack_max_algae, icepack_max_aero, icepack_max_fe
      18             :       use icepack_intfc, only: icepack_max_nbtrcr
      19             :       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
      20             :       use icepack_intfc, only: icepack_init_tracer_sizes, icepack_init_tracer_flags
      21             :       use icepack_intfc, only: icepack_init_tracer_indices
      22             :       use icepack_intfc, only: icepack_init_parameters
      23             :       use icepack_intfc, only: icepack_query_tracer_sizes, icepack_query_tracer_flags
      24             :       use icepack_intfc, only: icepack_query_tracer_indices
      25             :       use icepack_intfc, only: icepack_query_parameters
      26             :       use icepack_intfc, only: icepack_init_zbgc
      27             :       use icepack_intfc, only: icepack_init_thermo
      28             :       use icepack_intfc, only: icepack_step_radiation, icepack_init_orbit
      29             :       use icepack_intfc, only: icepack_init_bgc, icepack_init_zsalinity
      30             :       use icepack_intfc, only: icepack_init_ocean_bio, icepack_load_ocean_bio_array
      31             :       use icepack_intfc, only: icepack_init_hbrine
      32             :       use icedrv_system, only: icedrv_system_abort
      33             : 
      34             :       implicit none
      35             : 
      36             :       private
      37             :       public :: init_thermo_vertical, init_shortwave, &
      38             :                 init_bgc, init_hbrine, init_zbgc
      39             : 
      40             : !=======================================================================
      41             : 
      42             :       contains
      43             : 
      44             : !=======================================================================
      45             : !
      46             : ! Initialize the vertical profile of ice salinity and melting temperature.
      47             : !
      48             : ! authors: C. M. Bitz, UW
      49             : !          William H. Lipscomb, LANL
      50             : 
      51          46 :       subroutine init_thermo_vertical
      52             : 
      53             :       use icedrv_flux, only: salinz, Tmltz
      54             : 
      55             :       integer (kind=int_kind) :: &
      56             :          i,          &  ! horizontal indices
      57             :          k              ! ice layer index
      58             : 
      59             :       real (kind=dbl_kind), dimension(nilyr+1) :: &
      60         141 :          sprofile                         ! vertical salinity profile
      61             : 
      62             :       real (kind=dbl_kind) :: &
      63          17 :          depressT
      64             : 
      65             :       character(len=*), parameter :: subname='(init_thermo_vertical)'
      66             : 
      67             :       !-----------------------------------------------------------------
      68             :       ! query Icepack values
      69             :       !-----------------------------------------------------------------
      70             : 
      71          46 :       call icepack_query_parameters(depressT_out=depressT)
      72          46 :       call icepack_init_thermo(nilyr=nilyr, sprofile=sprofile)
      73          46 :       call icepack_warnings_flush(nu_diag)
      74          46 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
      75           0 :           file=__FILE__, line=__LINE__)
      76             : 
      77             :       !-----------------------------------------------------------------
      78             :       ! Prescibe vertical profile of salinity and melting temperature.
      79             :       ! Note this profile is only used for BL99 thermodynamics.
      80             :       !-----------------------------------------------------------------
      81             : 
      82         230 :       do i = 1, nx
      83        1558 :          do k = 1, nilyr+1
      84        1328 :             salinz(i,k) = sprofile(k)
      85        1512 :             Tmltz (i,k) = -salinz(i,k)*depressT
      86             :          enddo ! k
      87             :       enddo    ! i
      88             : 
      89          46 :       end subroutine init_thermo_vertical
      90             : 
      91             : !=======================================================================
      92             : !
      93             : !  Initialize shortwave
      94             : 
      95          46 :       subroutine init_shortwave
      96             : 
      97             :       use icedrv_arrays_column, only: fswpenln, Iswabsn, Sswabsn, albicen
      98             :       use icedrv_arrays_column, only: albsnon, alvdrn, alidrn, alvdfn, alidfn
      99             :       use icedrv_arrays_column, only: fswsfcn, ffracn, snowfracn
     100             :       use icedrv_arrays_column, only: fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf
     101             :       use icedrv_arrays_column, only: fswintn, albpndn, apeffn, trcrn_sw, dhsn
     102             :       use icedrv_arrays_column, only: kaer_tab, waer_tab, gaer_tab
     103             :       use icedrv_arrays_column, only: kaer_bc_tab, waer_bc_tab, gaer_bc_tab
     104             :       use icedrv_arrays_column, only: swgrid, igrid, bcenh
     105             :       use icedrv_calendar, only: istep1, dt, calendar_type
     106             :       use icedrv_calendar, only:    days_per_year, nextsw_cday, yday, sec
     107             :       use icedrv_system, only: icedrv_system_abort
     108             :       use icedrv_flux, only: alvdf, alidf, alvdr, alidr
     109             :       use icedrv_flux, only: alvdr_ai, alidr_ai, alvdf_ai, alidf_ai
     110             :       use icedrv_flux, only: swvdr, swvdf, swidr, swidf, scale_factor, snowfrac
     111             :       use icedrv_flux, only: albice, albsno, albpnd, apeff_ai, coszen, fsnow
     112             :       use icedrv_init, only: tlat, tlon, tmask
     113             :       use icedrv_restart_shared, only: restart
     114             :       use icedrv_state, only: aicen, vicen, vsnon, trcrn
     115             : 
     116             :       integer (kind=int_kind) :: &
     117             :          i, k         , & ! horizontal indices
     118             :          n                ! thickness category index
     119             : 
     120             :       real (kind=dbl_kind) :: &
     121          17 :          netsw            ! flag for shortwave radiation presence
     122             : 
     123             :       logical (kind=log_kind) :: &
     124             :          l_print_point, & ! flag to print designated grid point diagnostics
     125             :          dEdd_algae,    & ! from icepack
     126             :          modal_aero       ! from icepack
     127             : 
     128             :       character (len=char_len) :: &
     129             :          shortwave        ! from icepack
     130             : 
     131             :       real (kind=dbl_kind), dimension(ncat) :: &
     132          98 :          fbri             ! brine height to ice thickness
     133             : 
     134             :       real (kind=dbl_kind), allocatable, dimension(:,:) :: &
     135          46 :          ztrcr_sw
     136             : 
     137             :       logical (kind=log_kind) :: tr_brine, tr_zaero, tr_bgc_N
     138             :       integer (kind=int_kind) :: nt_alvl, nt_apnd, nt_hpnd, nt_ipnd, nt_aero, &
     139             :          nt_fbri, nt_tsfc, ntrcr, nbtrcr_sw, nlt_chl_sw
     140             :       integer (kind=int_kind), dimension(icepack_max_aero) :: nlt_zaero_sw
     141             :       integer (kind=int_kind), dimension(icepack_max_aero) :: nt_zaero
     142             :       integer (kind=int_kind), dimension(icepack_max_algae) :: nt_bgc_N
     143          17 :       real (kind=dbl_kind) :: puny
     144             : 
     145             :       character(len=*), parameter :: subname='(init_shortwave)'
     146             : 
     147             :       !-----------------------------------------------------------------
     148             :       ! query Icepack values
     149             :       !-----------------------------------------------------------------
     150             : 
     151          46 :          call icepack_query_parameters(puny_out=puny)
     152          46 :          call icepack_query_parameters(shortwave_out=shortwave)
     153          46 :          call icepack_query_parameters(dEdd_algae_out=dEdd_algae)
     154          46 :          call icepack_query_parameters(modal_aero_out=modal_aero)
     155             :          call icepack_query_tracer_sizes(ntrcr_out=ntrcr, &
     156          46 :               nbtrcr_sw_out=nbtrcr_sw)
     157             :          call icepack_query_tracer_flags(tr_brine_out=tr_brine, &
     158          46 :               tr_zaero_out=tr_zaero, tr_bgc_N_out=tr_bgc_N)
     159             :          call icepack_query_tracer_indices(nt_alvl_out=nt_alvl, &
     160             :               nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, &
     161             :               nt_ipnd_out=nt_ipnd, nt_aero_out=nt_aero, &
     162             :               nt_fbri_out=nt_fbri, nt_tsfc_out=nt_tsfc, &
     163             :               nt_bgc_N_out=nt_bgc_N, nt_zaero_out=nt_zaero, &
     164          46 :               nlt_chl_sw_out=nlt_chl_sw, nlt_zaero_sw_out=nlt_zaero_sw)
     165          46 :          call icepack_warnings_flush(nu_diag)
     166          46 :          if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     167           0 :              file=__FILE__,line= __LINE__)
     168             : 
     169             :       !-----------------------------------------------------------------
     170             :       ! Initialize
     171             :       !-----------------------------------------------------------------
     172             : 
     173          46 :          allocate(ztrcr_sw(nbtrcr_sw, ncat))
     174             : 
     175          46 :          fswpenln(:,:,:) = c0
     176          46 :          Iswabsn(:,:,:) = c0
     177          46 :          Sswabsn(:,:,:) = c0
     178             : 
     179         230 :          do i = 1, nx
     180             : 
     181         184 :             l_print_point = .false.
     182             : 
     183         184 :             alvdf(i) = c0
     184         184 :             alidf(i) = c0
     185         184 :             alvdr(i) = c0
     186         184 :             alidr(i) = c0
     187         184 :             alvdr_ai(i) = c0
     188         184 :             alidr_ai(i) = c0
     189         184 :             alvdf_ai(i) = c0
     190         184 :             alidf_ai(i) = c0
     191         184 :             albice(i) = c0
     192         184 :             albsno(i) = c0
     193         184 :             albpnd(i) = c0
     194         184 :             snowfrac(i) = c0
     195         184 :             apeff_ai(i) = c0
     196             : 
     197        1102 :             do n = 1, ncat
     198         872 :                alvdrn(i,n) = c0
     199         872 :                alidrn(i,n) = c0
     200         872 :                alvdfn(i,n) = c0
     201         872 :                alidfn(i,n) = c0
     202         872 :                fswsfcn(i,n) = c0
     203         872 :                fswintn(i,n) = c0
     204         872 :                fswthrun(i,n) = c0
     205         872 :                fswthrun_vdr(i,n) = c0
     206         872 :                fswthrun_vdf(i,n) = c0
     207         872 :                fswthrun_idr(i,n) = c0
     208        1056 :                fswthrun_idf(i,n) = c0
     209             :             enddo   ! ncat
     210             : 
     211             :          enddo
     212             : 
     213         230 :          do i = 1, nx
     214             : 
     215         184 :             if (trim(shortwave) == 'dEdd') then ! delta Eddington
     216             : 
     217             :                ! initialize orbital parameters
     218             :                ! These come from the driver in the coupled model.
     219         160 :                call icepack_warnings_flush(nu_diag)
     220         160 :                call icepack_init_orbit()
     221         160 :                call icepack_warnings_flush(nu_diag)
     222         160 :                if (icepack_warnings_aborted()) &
     223           0 :                   call icedrv_system_abort(i, istep1, subname, __FILE__, __LINE__)
     224             :             endif
     225             : 
     226         184 :             fbri(:) = c0
     227        1056 :             ztrcr_sw(:,:) = c0
     228        1056 :             do n = 1, ncat
     229        1056 :                if (tr_brine)  fbri(n) = trcrn(i,nt_fbri,n)
     230             :             enddo
     231             : 
     232         184 :             if (tmask(i)) then
     233             :             call icepack_step_radiation (              &
     234             :                          dt=dt,           ncat=ncat,           &
     235             :                          nblyr=nblyr,                          &
     236             :                          nilyr=nilyr,     nslyr=nslyr,         &
     237             :                          dEdd_algae=dEdd_algae,                &
     238             :                          swgrid=swgrid(:),                     &
     239             :                          igrid=igrid(:),                       &
     240             :                          fbri=fbri(:),                         &
     241             :                          aicen=aicen(i,:),                     &
     242             :                          vicen=vicen(i,:),                     &
     243             :                          vsnon=vsnon(i,:),                     &
     244             :                          Tsfcn=trcrn(i,nt_Tsfc,:),             &
     245             :                          alvln=trcrn(i,nt_alvl,:),             &
     246             :                          apndn=trcrn(i,nt_apnd,:),             &
     247             :                          hpndn=trcrn(i,nt_hpnd,:),             &
     248             :                          ipndn=trcrn(i,nt_ipnd,:),             &
     249           0 :                          aeron=trcrn(i,nt_aero:nt_aero+4*n_aero-1,:), &
     250           0 :                          bgcNn=trcrn(i,nt_bgc_N(1):nt_bgc_N(1)+n_algae*(nblyr+3)-1,:), &
     251           0 :                          zaeron=trcrn(i,nt_zaero(1):nt_zaero(1)+n_zaero*(nblyr+3)-1,:), &
     252             :                          trcrn_bgcsw=ztrcr_sw,                 &
     253          51 :                          TLAT=TLAT(i), TLON=TLON(i),           &
     254             :                          calendar_type=calendar_type,          &
     255             :                          days_per_year=days_per_year,          &
     256             :                          nextsw_cday=nextsw_cday, yday=yday, sec=sec,      &
     257             :                          kaer_tab=kaer_tab, kaer_bc_tab=kaer_bc_tab(:,:),  &
     258             :                          waer_tab=waer_tab, waer_bc_tab=waer_bc_tab(:,:),  &
     259             :                          gaer_tab=gaer_tab, gaer_bc_tab=gaer_bc_tab(:,:),  &
     260             :                          bcenh=bcenh(:,:,:),                               &
     261             :                          modal_aero=modal_aero,                            &
     262          51 :                          swvdr=swvdr(i),         swvdf=swvdf(i),           &
     263          51 :                          swidr=swidr(i),         swidf=swidf(i),           &
     264          51 :                          coszen=coszen(i),       fsnow=fsnow(i),           &
     265             :                          alvdrn=alvdrn(i,:),     alvdfn=alvdfn(i,:),       &
     266             :                          alidrn=alidrn(i,:),     alidfn=alidfn(i,:),       &
     267             :                          fswsfcn=fswsfcn(i,:),   fswintn=fswintn(i,:),     &
     268             :                          fswthrun=fswthrun(i,:),                           &
     269             :                          fswthrun_vdr=fswthrun_vdr(i,:),                   &
     270             :                          fswthrun_vdf=fswthrun_vdf(i,:),                   &
     271             :                          fswthrun_idr=fswthrun_idr(i,:),                   &
     272             :                          fswthrun_idf=fswthrun_idf(i,:),                   &
     273             :                          fswpenln=fswpenln(i,:,:),                         &
     274             :                          Sswabsn=Sswabsn(i,:,:), Iswabsn=Iswabsn(i,:,:),   &
     275             :                          albicen=albicen(i,:),   albsnon=albsnon(i,:),     &
     276             :                          albpndn=albpndn(i,:),   apeffn=apeffn(i,:),       &
     277             :                          snowfracn=snowfracn(i,:),                         &
     278             :                          dhsn=dhsn(i,:),         ffracn=ffracn(i,:),       &
     279             :                          l_print_point=l_print_point,                      &
     280         240 :                          initonly = .true.)
     281             :             endif
     282         184 :             call icepack_warnings_flush(nu_diag)
     283         184 :             if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     284           0 :                file=__FILE__, line=__LINE__)
     285             :          
     286             :       !-----------------------------------------------------------------
     287             :       ! Define aerosol tracer on shortwave grid
     288             :       !-----------------------------------------------------------------
     289             : 
     290         230 :             if (dEdd_algae .and. (tr_zaero .or. tr_bgc_N)) then
     291           0 :                do n = 1, ncat
     292         184 :                   do k = 1, nbtrcr_sw
     293           0 :                      trcrn_sw(i,k,n) = ztrcr_sw(k,n)
     294             :                   enddo
     295             :                enddo
     296             :             endif
     297             : 
     298             :          enddo
     299             : 
     300             :       !-----------------------------------------------------------------
     301             :       ! Aggregate albedos 
     302             :       !-----------------------------------------------------------------
     303             : 
     304         264 :          do n = 1, ncat
     305        1136 :             do i = 1, nx
     306             :                
     307        1090 :                if (aicen(i,n) > puny) then
     308             :                   
     309         372 :                   alvdf(i) = alvdf(i) + alvdfn(i,n)*aicen(i,n)
     310         372 :                   alidf(i) = alidf(i) + alidfn(i,n)*aicen(i,n)
     311         372 :                   alvdr(i) = alvdr(i) + alvdrn(i,n)*aicen(i,n)
     312         372 :                   alidr(i) = alidr(i) + alidrn(i,n)*aicen(i,n)
     313             :                   
     314         372 :                   netsw = swvdr(i) + swidr(i) + swvdf(i) + swidf(i)
     315         372 :                   if (netsw > puny) then ! sun above horizon
     316         227 :                      albice(i) = albice(i) + albicen(i,n)*aicen(i,n)
     317         227 :                      albsno(i) = albsno(i) + albsnon(i,n)*aicen(i,n)
     318         227 :                      albpnd(i) = albpnd(i) + albpndn(i,n)*aicen(i,n)
     319             :                   endif
     320             :                   
     321         372 :                   apeff_ai(i) = apeff_ai(i) + apeffn(i,n)*aicen(i,n)
     322         372 :                   snowfrac(i) = snowfrac(i) + snowfracn(i,n)*aicen(i,n)
     323             :                
     324             :                endif ! aicen > puny
     325             :             enddo  ! i
     326             :          enddo   ! ncat
     327             : 
     328         230 :          do i = 1, nx
     329             : 
     330             :       !----------------------------------------------------------------
     331             :       ! Store grid box mean albedos and fluxes before scaling by aice
     332             :       !----------------------------------------------------------------
     333             : 
     334         184 :             alvdf_ai  (i) = alvdf  (i)
     335         184 :             alidf_ai  (i) = alidf  (i)
     336         184 :             alvdr_ai  (i) = alvdr  (i)
     337         184 :             alidr_ai  (i) = alidr  (i)
     338             :             
     339             :       !----------------------------------------------------------------
     340             :       ! Save net shortwave for scaling factor in scale_factor
     341             :       !----------------------------------------------------------------
     342         230 :             if (.not. restart) then
     343          64 :                scale_factor(i) = &
     344          64 :                        swvdr(i)*(c1 - alvdr_ai(i)) &
     345          64 :                      + swvdf(i)*(c1 - alvdf_ai(i)) &
     346          64 :                      + swidr(i)*(c1 - alidr_ai(i)) &
     347         124 :                      + swidf(i)*(c1 - alidf_ai(i))
     348             :             endif
     349             : 
     350             :          enddo ! i
     351             : 
     352          46 :          deallocate(ztrcr_sw)
     353             : 
     354          46 :       end subroutine init_shortwave
     355             : 
     356             : !=======================================================================
     357             : 
     358             : !  Initialize vertical profile for biogeochemistry
     359             : 
     360           5 :       subroutine init_bgc() 
     361             : 
     362             :       use icedrv_arrays_column, only: zfswin, trcrn_sw
     363             :       use icedrv_arrays_column, only: ocean_bio_all, ice_bio_net, snow_bio_net
     364             :       use icedrv_arrays_column, only: cgrid, igrid, bphi, iDi, bTiz, iki
     365             :       use icedrv_arrays_column, only: Rayleigh_criteria, Rayleigh_real
     366             :       use icedrv_calendar,  only: istep1
     367             :       use icedrv_system, only: icedrv_system_abort
     368             :       use icedrv_flux, only: sss, nit, amm, sil, dmsp, dms, algalN, &
     369             :           doc, don, dic, fed, fep, zaeros, hum
     370             :       use icedrv_forcing_bgc, only:  get_forcing_bgc
     371             :       use icedrv_state, only: trcrn
     372             : 
     373             :       ! local variables
     374             : 
     375             :       integer (kind=int_kind) :: &
     376             :          i                , & ! horizontal indices
     377             :          k                , & ! vertical index
     378             :          n                    ! category index
     379             : 
     380             :       integer (kind=int_kind) :: &
     381             :          max_nbtrcr, max_algae, max_don, max_doc, max_dic, max_aero, max_fe
     382             : 
     383             :       logical (kind=log_kind) :: &
     384             :          RayleighC, &
     385             :          solve_zsal
     386             : 
     387             :       real(kind=dbl_kind) :: &
     388           3 :          RayleighR
     389             : 
     390             :       real(kind=dbl_kind), dimension(max_ntrcr,ncat) :: &
     391        2923 :          trcrn_bgc 
     392             :       
     393             :       real(kind=dbl_kind), dimension(nilyr,ncat) :: &
     394         123 :          sicen    
     395             : 
     396             :       integer (kind=int_kind) :: &
     397             :          nbtrcr, ntrcr, ntrcr_o, &
     398             :          nt_sice, nt_bgc_S
     399             : 
     400             :       character(len=*), parameter :: subname='(init_bgc)'
     401             : 
     402             :       ! Initialize
     403             : 
     404           5 :       bphi(:,:,:) = c0   ! initial porosity for no ice 
     405           5 :       iDi (:,:,:) = c0   ! interface diffusivity
     406           5 :       bTiz(:,:,:) = c0   ! initial bio grid ice temperature
     407           5 :       iki (:,:,:) = c0   ! permeability
     408             : 
     409           5 :       ocean_bio_all(:,:)   = c0
     410           5 :       ice_bio_net  (:,:)   = c0 ! integrated ice tracer conc (mmol/m^2 or mg/m^2) 
     411           5 :       snow_bio_net (:,:)   = c0 ! integrated snow tracer conc (mmol/m^2 or mg/m^2)
     412           5 :       zfswin       (:,:,:) = c0 ! shortwave flux on bio grid
     413           5 :       trcrn_sw     (:,:,:) = c0 ! tracers active in the shortwave calculation
     414           5 :       trcrn_bgc    (:,:)   = c0
     415             : 
     416             :       !-----------------------------------------------------------------
     417             :       ! query Icepack values
     418             :       !-----------------------------------------------------------------
     419             : 
     420           5 :       call icepack_query_parameters(solve_zsal_out=solve_zsal)
     421             :       call icepack_query_tracer_sizes(max_nbtrcr_out=max_nbtrcr, &
     422             :            max_algae_out=max_algae, max_don_out=max_don, max_doc_out=max_doc, &
     423           5 :            max_dic_out=max_dic, max_aero_out=max_aero, max_fe_out=max_fe)
     424             :       call icepack_query_tracer_sizes(nbtrcr_out=nbtrcr, ntrcr_out=ntrcr, &
     425           5 :            ntrcr_o_out=ntrcr_o)
     426           5 :       call icepack_query_tracer_indices(nt_sice_out=nt_sice, nt_bgc_S_out=nt_bgc_S)
     427           5 :       call icepack_warnings_flush(nu_diag)
     428           5 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     429           0 :           file=__FILE__,line= __LINE__)
     430             : 
     431             :       !-----------------------------------------------------------------
     432             :       ! zsalinity initialization
     433             :       !-----------------------------------------------------------------
     434             :       
     435           5 :       if (solve_zsal) then
     436           0 :          do i = 1, nx
     437             :             call icepack_init_zsalinity(ncat=ncat, nblyr=nblyr,        &
     438             :                                         ntrcr_o           = ntrcr_o,   &
     439             :                                         Rayleigh_criteria = RayleighC, &
     440             :                                         Rayleigh_real     = RayleighR, &
     441             :                                         trcrn_bgc         = trcrn_bgc, &
     442             :                                         nt_bgc_S          = nt_bgc_S,  &
     443           0 :                                         sss               = sss(i))
     444           0 :             Rayleigh_real    (i) = RayleighR
     445           0 :             Rayleigh_criteria(i) = RayleighC
     446           0 :             do n = 1,ncat
     447           0 :                do k  = 1, nblyr
     448           0 :                   trcrn(i,nt_bgc_S+k-1,n) = trcrn_bgc(nt_bgc_S-1+k-ntrcr_o,n)
     449             :                enddo
     450             :             enddo
     451             :          enddo      ! i
     452           0 :          call icepack_warnings_flush(nu_diag)
     453           0 :          if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     454           0 :              file=__FILE__, line=__LINE__)
     455             :       endif ! solve_zsal
     456             : 
     457             :       !-----------------------------------------------------------------
     458             :       ! biogeochemistry initialization
     459             :       !-----------------------------------------------------------------
     460             : 
     461             :       !-----------------------------------------------------------------
     462             :       ! Initial Ocean Values if not coupled to the ocean bgc
     463             :       !-----------------------------------------------------------------
     464          25 :       do i = 1, nx
     465             :          call icepack_init_ocean_bio ( &
     466          12 :                       amm=amm(i),   dmsp=dmsp(i),  dms=dms(i),   &
     467             :                       doc=doc(i,:), dic =dic(i,:), don=don(i,:), &
     468          12 :                       fed=fed(i,:), fep =fep(i,:), hum=hum(i),   &
     469          12 :                       nit=nit(i),   sil =sil(i),                 &
     470             :                       zaeros=zaeros(i,:),                        &
     471             :                       algalN=algalN(i,:),                        &
     472             :                       max_dic=max_dic, max_don =max_don,         &
     473          37 :                       max_fe =max_fe,  max_aero=max_aero)
     474             :       enddo  ! i
     475           5 :       call icepack_warnings_flush(nu_diag)
     476           5 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     477           0 :           file=__FILE__, line=__LINE__)
     478             : 
     479           5 :       call get_forcing_bgc                          ! defines nit and sil
     480             : 
     481          25 :       do i = 1, nx
     482             : 
     483         120 :          do n = 1, ncat
     484         800 :             do k = 1, nilyr
     485         800 :                sicen(k,n) = trcrn(i,nt_sice+k-1,n)
     486             :             enddo
     487       17260 :             do k = ntrcr_o+1, ntrcr
     488       17240 :                trcrn_bgc(k-ntrcr_o,n) = trcrn(i,k,n)
     489             :             enddo
     490             :          enddo
     491             : 
     492             :          call icepack_load_ocean_bio_array(max_nbtrcr=max_nbtrcr,             &
     493             :                       max_algae=max_algae, max_don=max_don,  max_doc=max_doc, &
     494             :                       max_aero =max_aero,  max_dic=max_dic,  max_fe =max_fe,  &
     495          12 :                       nit =nit(i),   amm=amm(i),   sil   =sil(i),             &
     496          12 :                       dmsp=dmsp(i),  dms=dms(i),   algalN=algalN(i,:),        &
     497             :                       doc =doc(i,:), don=don(i,:), dic   =dic(i,:),           &  
     498             :                       fed =fed(i,:), fep=fep(i,:), zaeros=zaeros(i,:),        &
     499          20 :                       ocean_bio_all=ocean_bio_all(i,:),  hum=hum(i))
     500          20 :          call icepack_warnings_flush(nu_diag)
     501          20 :          if (icepack_warnings_aborted()) call icedrv_system_abort(i, istep1, subname, &
     502           5 :              __FILE__, __LINE__)
     503             : 
     504             :       enddo  ! i
     505             : 
     506          25 :       do i = 1, nx
     507             :          call icepack_init_bgc(ncat=ncat, nblyr=nblyr, nilyr=nilyr, ntrcr_o=ntrcr_o, &
     508             :                       cgrid=cgrid, igrid=igrid, ntrcr=ntrcr, nbtrcr=nbtrcr,          &
     509             :                       sicen=sicen(:,:), trcrn=trcrn_bgc(:,:),                        &
     510          20 :                       sss=sss(i), ocean_bio_all=ocean_bio_all(i,:))
     511          20 :          call icepack_warnings_flush(nu_diag)
     512          20 :          if (icepack_warnings_aborted()) call icedrv_system_abort(i, istep1, subname, &
     513           5 :              __FILE__, __LINE__)
     514             :       enddo  ! i
     515             : 
     516           5 :       end subroutine init_bgc
     517             : 
     518             : !=======================================================================
     519             : 
     520             : !  Initialize brine height tracer
     521             : 
     522           5 :       subroutine init_hbrine()
     523             : 
     524             :       use icedrv_arrays_column, only: first_ice, bgrid, igrid, cgrid
     525             :       use icedrv_arrays_column, only: icgrid, swgrid
     526             :       use icedrv_state, only: trcrn
     527             : 
     528           3 :       real (kind=dbl_kind) :: phi_snow
     529             :       integer (kind=int_kind) :: nt_fbri
     530             :       logical (kind=log_kind) :: tr_brine
     531             :       character(len=*), parameter :: subname='(init_hbrine)'
     532             : 
     533             :       !-----------------------------------------------------------------
     534             :       ! query Icepack values
     535             :       !-----------------------------------------------------------------
     536             : 
     537           5 :       call icepack_query_parameters(phi_snow_out=phi_snow)
     538           5 :       call icepack_query_tracer_flags(tr_brine_out=tr_brine)
     539           5 :       call icepack_query_tracer_indices(nt_fbri_out=nt_fbri)
     540           5 :       call icepack_warnings_flush(nu_diag)
     541           5 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     542           0 :           file=__FILE__,line= __LINE__)
     543             : 
     544             :       !-----------------------------------------------------------------
     545             : 
     546             :       call icepack_init_hbrine(bgrid=bgrid, igrid=igrid, cgrid=cgrid, icgrid=icgrid, &
     547           5 :            swgrid=swgrid, nblyr=nblyr, nilyr=nilyr, phi_snow=phi_snow)
     548           5 :       call icepack_warnings_flush(nu_diag)
     549           5 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     550           0 :           file=__FILE__, line=__LINE__)
     551             : 
     552           5 :       call icepack_init_parameters(phi_snow_in=phi_snow)
     553           5 :       call icepack_warnings_flush(nu_diag)
     554           5 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     555           0 :           file=__FILE__,line= __LINE__)
     556             : 
     557         130 :       first_ice(:,:) = .true.            
     558         130 :       if (tr_brine) trcrn(:,nt_fbri,:) = c1
     559             : 
     560           5 :       end subroutine init_hbrine
     561             : 
     562             : !=======================================================================
     563             : 
     564             : ! Namelist variables, set to default values; may be altered at run time
     565             : ! 
     566             : ! author Elizabeth C. Hunke, LANL
     567             : !        Nicole Jeffery, LANL
     568             : 
     569          46 :       subroutine init_zbgc
     570             : 
     571             :       use icedrv_state, only: trcr_base, trcr_depend, n_trcr_strata
     572             :       use icedrv_state, only: nt_strata
     573             :       use icedrv_forcing, only: bgc_data_type
     574             : 
     575             :       character (len=char_len) :: &
     576             :          shortwave        ! from icepack
     577             : 
     578             :       integer (kind=int_kind) :: &
     579             :          ntrcr,         nbtrcr,       nbtrcr_sw,    &
     580             :          ntrcr_o,       nt_fbri,      &
     581             :          nt_bgc_Nit,    nt_bgc_Am,    nt_bgc_Sil,   &
     582             :          nt_bgc_DMS,    nt_bgc_PON,   nt_bgc_S,     &
     583             :          nt_bgc_DMSPp,  nt_bgc_DMSPd, &
     584             :          nt_zbgc_frac,  nlt_chl_sw,   &
     585             :          nlt_bgc_Nit,   nlt_bgc_Am,   nlt_bgc_Sil,  &
     586             :          nlt_bgc_DMS,   nlt_bgc_DMSPp,nlt_bgc_DMSPd,&
     587             :          nlt_bgc_PON,   nt_bgc_hum,   nlt_bgc_hum
     588             : 
     589             :       integer (kind=int_kind), dimension(icepack_max_aero) :: &
     590             :          nlt_zaero_sw       ! points to aerosol in trcrn_sw
     591             : 
     592             :       integer (kind=int_kind), dimension(icepack_max_algae) :: &
     593             :          nlt_bgc_N      , & ! algae
     594             :          nlt_bgc_C      , & !
     595             :          nlt_bgc_chl
     596             : 
     597             :       integer (kind=int_kind), dimension(icepack_max_doc) :: &
     598             :          nlt_bgc_DOC        ! disolved organic carbon
     599             : 
     600             :       integer (kind=int_kind), dimension(icepack_max_don) :: &
     601             :          nlt_bgc_DON        !
     602             : 
     603             :       integer (kind=int_kind), dimension(icepack_max_dic) :: &
     604             :          nlt_bgc_DIC        ! disolved inorganic carbon
     605             : 
     606             :       integer (kind=int_kind), dimension(icepack_max_fe) :: &
     607             :          nlt_bgc_Fed    , & !
     608             :          nlt_bgc_Fep        !
     609             : 
     610             :       integer (kind=int_kind), dimension(icepack_max_aero) :: &
     611             :          nlt_zaero          ! non-reacting layer aerosols
     612             : 
     613             :       integer (kind=int_kind), dimension(icepack_max_algae) :: &
     614             :          nt_bgc_N , & ! diatoms, phaeocystis, pico/small
     615             :          nt_bgc_C , & ! diatoms, phaeocystis, pico/small
     616             :          nt_bgc_chl   ! diatoms, phaeocystis, pico/small
     617             : 
     618             :       integer (kind=int_kind), dimension(icepack_max_doc) :: &
     619             :          nt_bgc_DOC         !  dissolved organic carbon
     620             : 
     621             :       integer (kind=int_kind), dimension(icepack_max_don) :: &
     622             :          nt_bgc_DON         !  dissolved organic nitrogen
     623             : 
     624             :       integer (kind=int_kind), dimension(icepack_max_dic) :: &
     625             :          nt_bgc_DIC         !  dissolved inorganic carbon
     626             : 
     627             :       integer (kind=int_kind), dimension(icepack_max_fe) :: &
     628             :          nt_bgc_Fed     , & !  dissolved iron
     629             :          nt_bgc_Fep         !  particulate iron
     630             : 
     631             :       integer (kind=int_kind), dimension(icepack_max_aero) :: &
     632             :          nt_zaero       , & !  black carbon and other aerosols
     633             :          nt_zaero_sw        !  for shortwave calculation
     634             : 
     635             :       integer (kind=int_kind), dimension(icepack_max_nbtrcr) :: &
     636             :          bio_index_o        ! relates nlt_bgc_NO to ocean concentration index
     637             : 
     638             :       integer (kind=int_kind), dimension(icepack_max_nbtrcr) :: &
     639             :          bio_index          ! relates bio indices nlt_bgc_N to nt_bgc_N
     640             : 
     641             :       logical (kind=log_kind) :: &
     642             :           tr_brine, &
     643             :           tr_bgc_Nit,    tr_bgc_Am,    tr_bgc_Sil,   &
     644             :           tr_bgc_DMS,    tr_bgc_PON,                 &
     645             :           tr_bgc_N,      tr_bgc_C,     tr_bgc_chl,   &
     646             :           tr_bgc_DON,    tr_bgc_Fe,    tr_zaero,     &
     647             :           tr_bgc_hum,    tr_aero
     648             :  
     649             :       integer (kind=int_kind) :: &
     650             :           ktherm
     651             : 
     652             :       logical (kind=log_kind) :: &
     653             :           solve_zsal, skl_bgc, z_tracers, scale_bgc, solve_zbgc, dEdd_algae, &
     654             :           modal_aero, restore_bgc
     655             : 
     656             :       character (char_len) :: &
     657             :           bgc_flux_type
     658             : 
     659             :       real (kind=dbl_kind) :: &
     660          51 :           grid_o, l_sk, initbio_frac, &
     661          51 :           frazil_scav, grid_oS, l_skS, &
     662          17 :           phi_snow, &
     663          17 :           ratio_Si2N_diatoms , ratio_Si2N_sp      , ratio_Si2N_phaeo   ,  &
     664          17 :           ratio_S2N_diatoms  , ratio_S2N_sp       , ratio_S2N_phaeo    ,  &
     665          17 :           ratio_Fe2C_diatoms , ratio_Fe2C_sp      , ratio_Fe2C_phaeo   ,  &
     666          17 :           ratio_Fe2N_diatoms , ratio_Fe2N_sp      , ratio_Fe2N_phaeo   ,  &
     667          17 :           ratio_Fe2DON       , ratio_Fe2DOC_s     , ratio_Fe2DOC_l     ,  &
     668          34 :           fr_resp            , tau_min            , tau_max            ,  &
     669          51 :           algal_vel          , R_dFe2dust         , dustFe_sol         ,  &
     670          17 :           chlabs_diatoms     , chlabs_sp          , chlabs_phaeo       ,  &
     671          17 :           alpha2max_low_diatoms,alpha2max_low_sp  , alpha2max_low_phaeo,  &
     672          17 :           beta2max_diatoms   , beta2max_sp        , beta2max_phaeo     ,  &
     673          17 :           mu_max_diatoms     , mu_max_sp          , mu_max_phaeo       ,  &
     674          17 :           grow_Tdep_diatoms  , grow_Tdep_sp       , grow_Tdep_phaeo    ,  &
     675          51 :           fr_graze_diatoms   , fr_graze_sp        , fr_graze_phaeo     ,  &
     676          17 :           mort_pre_diatoms   , mort_pre_sp        , mort_pre_phaeo     ,  &
     677          17 :           mort_Tdep_diatoms  , mort_Tdep_sp       , mort_Tdep_phaeo    ,  &
     678          17 :           k_exude_diatoms    , k_exude_sp         , k_exude_phaeo      ,  &
     679          17 :           K_Nit_diatoms      , K_Nit_sp           , K_Nit_phaeo        ,  &
     680          17 :           K_Am_diatoms       , K_Am_sp            , K_Am_phaeo         ,  &
     681          17 :           K_Sil_diatoms      , K_Sil_sp           , K_Sil_phaeo        ,  &
     682          17 :           K_Fe_diatoms       , K_Fe_sp            , K_Fe_phaeo         ,  &
     683          34 :           f_don_protein      , kn_bac_protein     , f_don_Am_protein   ,  &
     684          34 :           f_doc_s            , f_doc_l            , f_exude_s          ,  &
     685          34 :           f_exude_l          , k_bac_s            , k_bac_l            ,  &
     686          51 :           T_max              , fsal               , op_dep_min         ,  &
     687          51 :           fr_graze_s         , fr_graze_e         , fr_mort2min        ,  &
     688          51 :           fr_dFe             , k_nitrif           , t_iron_conv        ,  &
     689          34 :           max_loss           , max_dfe_doc1       , fr_resp_s          ,  &
     690          34 :           y_sk_DMS           , t_sk_conv          , t_sk_ox            ,  &
     691          17 :           algaltype_diatoms  , algaltype_sp       , algaltype_phaeo    ,  &
     692          51 :           nitratetype        , ammoniumtype       , silicatetype       ,  &
     693          34 :           dmspptype          , dmspdtype          , humtype            ,  &
     694          34 :           doctype_s          , doctype_l          , dontype_protein    ,  &
     695          51 :           fedtype_1          , feptype_1          , zaerotype_bc1      ,  &
     696          17 :           zaerotype_bc2      , zaerotype_dust1    , zaerotype_dust2    ,  &
     697          34 :           zaerotype_dust3    , zaerotype_dust4    , ratio_C2N_diatoms  ,  &
     698          34 :           ratio_C2N_sp       , ratio_C2N_phaeo    , ratio_chl2N_diatoms,  & 
     699          34 :           ratio_chl2N_sp     , ratio_chl2N_phaeo  , F_abs_chl_diatoms  ,  &
     700          34 :           F_abs_chl_sp       , F_abs_chl_phaeo    , ratio_C2N_proteins 
     701             : 
     702             :       real (kind=dbl_kind), dimension(icepack_max_dic) :: &
     703          34 :          dictype
     704             : 
     705             :       real (kind=dbl_kind), dimension(icepack_max_algae) :: &
     706          68 :          algaltype   ! tau_min for both retention and release
     707             : 
     708             :       real (kind=dbl_kind), dimension(icepack_max_doc) :: &
     709          68 :          doctype
     710             : 
     711             :       real (kind=dbl_kind), dimension(icepack_max_don) :: &
     712          34 :          dontype
     713             : 
     714             :       real (kind=dbl_kind), dimension(icepack_max_fe) :: &
     715          51 :          fedtype
     716             : 
     717             :       real (kind=dbl_kind), dimension(icepack_max_fe) :: &
     718          51 :          feptype
     719             : 
     720             :       real (kind=dbl_kind), dimension(icepack_max_aero) :: &
     721         119 :          zaerotype
     722             : 
     723             :       real (kind=dbl_kind), dimension(icepack_max_algae) :: &
     724          68 :          R_C2N     ,      & ! algal C to N (mole/mole)
     725          68 :          R_chl2N   ,      & ! 3 algal chlorophyll to N (mg/mmol)
     726          68 :          F_abs_chl          ! to scale absorption in Dedd
     727             : 
     728             :       real (kind=dbl_kind), dimension(icepack_max_don) :: &  ! increase compare to algal R_Fe2C
     729          34 :          R_C2N_DON
     730             : 
     731             :        real (kind=dbl_kind),  dimension(icepack_max_algae) :: &
     732          68 :          R_Si2N     , & ! algal Sil to N (mole/mole) 
     733          68 :          R_S2N      , & ! algal S to N (mole/mole)
     734             :          ! Marchetti et al 2006, 3 umol Fe/mol C for iron limited Pseudo-nitzschia
     735          68 :          R_Fe2C     , & ! algal Fe to carbon (umol/mmol)
     736          68 :          R_Fe2N         ! algal Fe to N (umol/mmol)
     737             : 
     738             :       real (kind=dbl_kind), dimension(icepack_max_don) :: & 
     739          34 :          R_Fe2DON       ! Fe to N of DON (nmol/umol)
     740             : 
     741             :       real (kind=dbl_kind), dimension(icepack_max_doc) :: &  
     742          68 :          R_Fe2DOC       ! Fe to C of DOC (nmol/umol)
     743             : 
     744             :       real (kind=dbl_kind), dimension(icepack_max_algae) :: &
     745          68 :          chlabs           , & ! chla absorption 1/m/(mg/m^3)
     746          68 :          alpha2max_low    , & ! light limitation (1/(W/m^2))
     747          68 :          beta2max         , & ! light inhibition (1/(W/m^2))
     748          68 :          mu_max           , & ! maximum growth rate (1/d)
     749          68 :          grow_Tdep        , & ! T dependence of growth (1/C)
     750          68 :          fr_graze         , & ! fraction of algae grazed
     751          68 :          mort_pre         , & ! mortality (1/day)
     752          68 :          mort_Tdep        , & ! T dependence of mortality (1/C)
     753          68 :          k_exude          , & ! algal carbon  exudation rate (1/d)
     754          68 :          K_Nit            , & ! nitrate half saturation (mmol/m^3) 
     755          68 :          K_Am             , & ! ammonium half saturation (mmol/m^3) 
     756          68 :          K_Sil            , & ! silicon half saturation (mmol/m^3)
     757          68 :          K_Fe                 ! iron half saturation  or micromol/m^3
     758             :             
     759             :       real (kind=dbl_kind), dimension(icepack_max_DON) :: &
     760          34 :          f_don            , & ! fraction of spilled grazing to DON
     761          34 :          kn_bac           , & ! Bacterial degredation of DON (1/d)
     762          34 :          f_don_Am             ! fraction of remineralized DON to Am
     763             : 
     764             :       real (kind=dbl_kind), dimension(icepack_max_DOC) :: &
     765          68 :          f_exude          , & ! fraction of exuded carbon to each DOC pool
     766          68 :          k_bac                ! Bacterial degredation of DOC (1/d)    
     767             : 
     768             :       real (kind=dbl_kind), dimension(icepack_max_nbtrcr) :: &
     769         510 :          zbgc_frac_init,&! initializes mobile fraction
     770         510 :          bgc_tracer_type ! described tracer in mobile or stationary phases      
     771             :                          ! < 0 is purely mobile (eg. nitrate)
     772             :                          ! > 0 has timescales for transitions between 
     773             :                          ! phases based on whether the ice is melting or growing
     774             : 
     775             :      real (kind=dbl_kind), dimension(icepack_max_nbtrcr) :: &
     776         510 :          zbgc_init_frac, &   ! fraction of ocean tracer  concentration in new ice
     777         510 :          tau_ret,        &   ! retention timescale  (s), mobile to stationary phase
     778         510 :          tau_rel             ! release timescale    (s), stationary to mobile phase
     779             : 
     780             :       character (32) :: &
     781             :          nml_filename = 'icepack_in' ! namelist input file name
     782             : 
     783             :       integer (kind=int_kind) :: &
     784             :         nml_error, & ! namelist i/o error flag
     785             :         k, mm    , & ! loop index
     786             :         ntd      , & ! for tracer dependency calculation
     787             :         nk       , & !
     788             :         nt_depend
     789             : 
     790             :       character(len=*), parameter :: subname='(init_zbgc)'
     791             : 
     792             :       !------------------------------------------------------------
     793             :       !        Tracers have mobile and  stationary phases. 
     794             :       ! ice growth allows for retention, ice melt facilitates mobility
     795             :       ! bgc_tracer_type defines the exchange timescales between these phases
     796             :       ! -1 : entirely in the mobile phase, no exchange  (this is the default)
     797             :       !  0 : retention time scale is tau_min, release time scale is tau_max
     798             :       !  1 : retention time scale is tau_max, release time scale is tau_min
     799             :       ! 0.5: retention time scale is tau_min, release time scale is tau_min
     800             :       !  2 : retention time scale is tau_max, release time scale is tau_max
     801             :       !------------------------------------------------------------
     802             : 
     803             :       !-----------------------------------------------------------------
     804             :       ! namelist variables
     805             :       !-----------------------------------------------------------------
     806             : 
     807             :       namelist /zbgc_nml/  &
     808             :         tr_brine, tr_zaero, modal_aero, skl_bgc, &
     809             :         z_tracers, dEdd_algae, solve_zbgc, bgc_flux_type, &
     810             :         restore_bgc, scale_bgc, solve_zsal, bgc_data_type, &
     811             :         tr_bgc_Nit, tr_bgc_C, tr_bgc_chl, tr_bgc_Am, tr_bgc_Sil, &
     812             :         tr_bgc_DMS, tr_bgc_PON, tr_bgc_hum, tr_bgc_DON, tr_bgc_Fe, &
     813             :         grid_o, l_sk, grid_oS, &
     814             :         l_skS, phi_snow,  initbio_frac, frazil_scav, &
     815             :         ratio_Si2N_diatoms , ratio_Si2N_sp      , ratio_Si2N_phaeo   ,  &
     816             :         ratio_S2N_diatoms  , ratio_S2N_sp       , ratio_S2N_phaeo    ,  &
     817             :         ratio_Fe2C_diatoms , ratio_Fe2C_sp      , ratio_Fe2C_phaeo   ,  &
     818             :         ratio_Fe2N_diatoms , ratio_Fe2N_sp      , ratio_Fe2N_phaeo   ,  &
     819             :         ratio_Fe2DON       , ratio_Fe2DOC_s     , ratio_Fe2DOC_l     ,  &
     820             :         fr_resp            , tau_min            , tau_max            ,  &
     821             :         algal_vel          , R_dFe2dust         , dustFe_sol         ,  &
     822             :         chlabs_diatoms     , chlabs_sp          , chlabs_phaeo       ,  &
     823             :         alpha2max_low_diatoms,alpha2max_low_sp  , alpha2max_low_phaeo,  &
     824             :         beta2max_diatoms   , beta2max_sp        , beta2max_phaeo     ,  &
     825             :         mu_max_diatoms     , mu_max_sp          , mu_max_phaeo       ,  &
     826             :         grow_Tdep_diatoms  , grow_Tdep_sp       , grow_Tdep_phaeo    ,  &
     827             :         fr_graze_diatoms   , fr_graze_sp        , fr_graze_phaeo     ,  &
     828             :         mort_pre_diatoms   , mort_pre_sp        , mort_pre_phaeo     ,  &
     829             :         mort_Tdep_diatoms  , mort_Tdep_sp       , mort_Tdep_phaeo    ,  &
     830             :         k_exude_diatoms    , k_exude_sp         , k_exude_phaeo      ,  &
     831             :         K_Nit_diatoms      , K_Nit_sp           , K_Nit_phaeo        ,  &
     832             :         K_Am_diatoms       , K_Am_sp            , K_Am_phaeo         ,  &
     833             :         K_Sil_diatoms      , K_Sil_sp           , K_Sil_phaeo        ,  &
     834             :         K_Fe_diatoms       , K_Fe_sp            , K_Fe_phaeo         ,  &
     835             :         f_don_protein      , kn_bac_protein     , f_don_Am_protein   ,  &
     836             :         f_doc_s            , f_doc_l            , f_exude_s          ,  &
     837             :         f_exude_l          , k_bac_s            , k_bac_l            ,  &
     838             :         T_max              , fsal               , op_dep_min         ,  &
     839             :         fr_graze_s         , fr_graze_e         , fr_mort2min        ,  &
     840             :         fr_dFe             , k_nitrif           , t_iron_conv        ,  &
     841             :         max_loss           , max_dfe_doc1       , fr_resp_s          ,  &
     842             :         y_sk_DMS           , t_sk_conv          , t_sk_ox            ,  &
     843             :         algaltype_diatoms  , algaltype_sp       , algaltype_phaeo    ,  &
     844             :         nitratetype        , ammoniumtype       , silicatetype       ,  &
     845             :         dmspptype          , dmspdtype          , humtype            ,  &
     846             :         doctype_s          , doctype_l          , dontype_protein    ,  &
     847             :         fedtype_1          , feptype_1          , zaerotype_bc1      ,  &
     848             :         zaerotype_bc2      , zaerotype_dust1    , zaerotype_dust2    ,  &
     849             :         zaerotype_dust3    , zaerotype_dust4    , ratio_C2N_diatoms  ,  &
     850             :         ratio_C2N_sp       , ratio_C2N_phaeo    , ratio_chl2N_diatoms,  & 
     851             :         ratio_chl2N_sp     , ratio_chl2N_phaeo  , F_abs_chl_diatoms  ,  &
     852             :         F_abs_chl_sp       , F_abs_chl_phaeo    , ratio_C2N_proteins 
     853             : 
     854             :       !-----------------------------------------------------------------
     855             :       ! query Icepack values
     856             :       !-----------------------------------------------------------------
     857             : 
     858          46 :       call icepack_query_tracer_sizes(ntrcr_out=ntrcr)
     859          46 :       call icepack_query_tracer_flags(tr_aero_out=tr_aero)
     860             :       call icepack_query_parameters(ktherm_out=ktherm, shortwave_out=shortwave, &
     861             :            scale_bgc_out=scale_bgc, skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, &
     862             :            dEdd_algae_out=dEdd_algae, solve_zbgc_out=solve_zbgc, phi_snow_out=phi_snow, &
     863             :            bgc_flux_type_out=bgc_flux_type, grid_o_out=grid_o, l_sk_out=l_sk, &
     864             :            initbio_frac_out=initbio_frac, frazil_scav_out=frazil_scav, &
     865             :            algal_vel_out=algal_vel, R_dFe2dust_out=R_dFe2dust, &
     866             :            dustFe_sol_out=dustFe_sol, T_max_out=T_max, fsal_out=fsal, &
     867             :            op_dep_min_out=op_dep_min, fr_graze_s_out=fr_graze_s, &
     868             :            fr_graze_e_out=fr_graze_e, fr_mort2min_out=fr_mort2min, &
     869             :            fr_dFe_out=fr_dFe, k_nitrif_out=k_nitrif, t_iron_conv_out=t_iron_conv, &
     870             :            max_loss_out=max_loss, max_dfe_doc1_out=max_dfe_doc1, &
     871             :            fr_resp_s_out=fr_resp_s, y_sk_DMS_out=y_sk_DMS, t_sk_conv_out=t_sk_conv, &
     872          46 :            t_sk_ox_out=t_sk_ox)
     873          46 :       call icepack_warnings_flush(nu_diag)
     874          46 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
     875           0 :           file=__FILE__, line=__LINE__)
     876             : 
     877             :       !-----------------------------------------------------------------
     878             :       ! default values
     879             :       !-----------------------------------------------------------------
     880          46 :       tr_brine        = .false.  ! brine height differs from ice height
     881          46 :       tr_zaero        = .false.  ! z aerosol tracers
     882          46 :       modal_aero      = .false.  ! use modal aerosol treatment of aerosols
     883          46 :       restore_bgc     = .false.  ! restore bgc if true
     884          46 :       solve_zsal      = .false.  ! update salinity tracer profile from solve_S_dt
     885          46 :       bgc_data_type   = 'default'! source of bgc data
     886          46 :       tr_bgc_PON      = .false.  !---------------------------------------------   
     887          46 :       tr_bgc_Nit      = .false.  ! biogeochemistry (skl or zbgc)
     888          46 :       tr_bgc_C        = .false.  ! if skl_bgc = .true. then skl
     889          46 :       tr_bgc_chl      = .false.  ! if z_tracers = .true. then vertically resolved
     890          46 :       tr_bgc_Sil      = .false.  ! if z_tracers + solve_zbgc = .true. then
     891          46 :       tr_bgc_Am       = .false.  ! vertically resolved with reactions  
     892          46 :       tr_bgc_DMS      = .false.  !------------------------------------------------
     893          46 :       tr_bgc_DON      = .false.  ! 
     894          46 :       tr_bgc_hum      = .false.  !
     895          46 :       tr_bgc_Fe       = .false.  ! 
     896          46 :       tr_bgc_N        = .true.   !
     897             : 
     898             :       ! z biology parameters  
     899          46 :       ratio_Si2N_diatoms = 1.8_dbl_kind    ! algal Si to N (mol/mol)                       
     900          46 :       ratio_Si2N_sp      = c0              ! diatoms, small plankton, phaeocystis
     901          46 :       ratio_Si2N_phaeo   = c0
     902          46 :       ratio_S2N_diatoms  = 0.03_dbl_kind   ! algal S  to N (mol/mol)
     903          46 :       ratio_S2N_sp       = 0.03_dbl_kind 
     904          46 :       ratio_S2N_phaeo    = 0.03_dbl_kind
     905          46 :       ratio_Fe2C_diatoms = 0.0033_dbl_kind ! algal Fe to C  (umol/mol)
     906          46 :       ratio_Fe2C_sp      = 0.0033_dbl_kind
     907          46 :       ratio_Fe2C_phaeo   = p1
     908          46 :       ratio_Fe2N_diatoms = 0.023_dbl_kind  ! algal Fe to N  (umol/mol)
     909          46 :       ratio_Fe2N_sp      = 0.023_dbl_kind
     910          46 :       ratio_Fe2N_phaeo   = 0.7_dbl_kind
     911          46 :       ratio_Fe2DON       = 0.023_dbl_kind  ! Fe to N of DON (nmol/umol)
     912          46 :       ratio_Fe2DOC_s     = p1              ! Fe to C of DOC (nmol/umol) saccharids
     913          46 :       ratio_Fe2DOC_l     = 0.033_dbl_kind  ! Fe to C of DOC (nmol/umol) lipids
     914          46 :       tau_min            = 5200.0_dbl_kind ! rapid mobile to stationary exchanges (s)
     915          46 :       tau_max            = 1.73e5_dbl_kind ! long time mobile to stationary exchanges (s)
     916          46 :       chlabs_diatoms     = 0.03_dbl_kind   ! chl absorption (1/m/(mg/m^3))
     917          46 :       chlabs_sp          = 0.01_dbl_kind
     918          46 :       chlabs_phaeo       = 0.05_dbl_kind
     919          46 :       alpha2max_low_diatoms = 0.8_dbl_kind ! light limitation (1/(W/m^2))  
     920          46 :       alpha2max_low_sp      = 0.67_dbl_kind
     921          46 :       alpha2max_low_phaeo   = 0.67_dbl_kind
     922          46 :       beta2max_diatoms   = 0.018_dbl_kind  ! light inhibition (1/(W/m^2))  
     923          46 :       beta2max_sp        = 0.0025_dbl_kind
     924          46 :       beta2max_phaeo     = 0.01_dbl_kind
     925          46 :       mu_max_diatoms     = 1.2_dbl_kind    ! maximum growth rate (1/day) 
     926          46 :       mu_max_sp          = 0.851_dbl_kind
     927          46 :       mu_max_phaeo       = 0.851_dbl_kind
     928          46 :       grow_Tdep_diatoms  = 0.06_dbl_kind ! Temperature dependence of growth (1/C)
     929          46 :       grow_Tdep_sp       = 0.06_dbl_kind
     930          46 :       grow_Tdep_phaeo    = 0.06_dbl_kind
     931          46 :       fr_graze_diatoms   = 0.01_dbl_kind ! Fraction grazed
     932          46 :       fr_graze_sp        = p1
     933          46 :       fr_graze_phaeo     = p1
     934          46 :       mort_pre_diatoms   = 0.007_dbl_kind! Mortality (1/day)
     935          46 :       mort_pre_sp        = 0.007_dbl_kind
     936          46 :       mort_pre_phaeo     = 0.007_dbl_kind
     937          46 :       mort_Tdep_diatoms  = 0.03_dbl_kind ! T dependence of mortality (1/C)
     938          46 :       mort_Tdep_sp       = 0.03_dbl_kind
     939          46 :       mort_Tdep_phaeo    = 0.03_dbl_kind
     940          46 :       k_exude_diatoms    = c0            ! algal exudation (1/d)
     941          46 :       k_exude_sp         = c0
     942          46 :       k_exude_phaeo      = c0
     943          46 :       K_Nit_diatoms      = c1            ! nitrate half saturation (mmol/m^3)
     944          46 :       K_Nit_sp           = c1
     945          46 :       K_Nit_phaeo        = c1
     946          46 :       K_Am_diatoms       = 0.3_dbl_kind  ! ammonium half saturation (mmol/m^3)
     947          46 :       K_Am_sp            = 0.3_dbl_kind
     948          46 :       K_Am_phaeo         = 0.3_dbl_kind
     949          46 :       K_Sil_diatoms      = 4.0_dbl_kind  ! silicate half saturation (mmol/m^3)
     950          46 :       K_Sil_sp           = c0
     951          46 :       K_Sil_phaeo        = c0
     952          46 :       K_Fe_diatoms       = c1            ! iron half saturation (nM)
     953          46 :       K_Fe_sp            = 0.2_dbl_kind
     954          46 :       K_Fe_phaeo         = p1
     955          46 :       f_don_protein      = 0.6_dbl_kind  ! fraction of spilled grazing to proteins           
     956          46 :       kn_bac_protein     = 0.03_dbl_kind ! Bacterial degredation of DON (1/d)                
     957          46 :       f_don_Am_protein   = 0.25_dbl_kind ! fraction of remineralized DON to ammonium         
     958          46 :       f_doc_s            = 0.4_dbl_kind  ! fraction of mortality to DOC 
     959          46 :       f_doc_l            = 0.4_dbl_kind
     960          46 :       f_exude_s          = c1            ! fraction of exudation to DOC
     961          46 :       f_exude_l          = c1
     962          46 :       k_bac_s            = 0.03_dbl_kind ! Bacterial degredation of DOC (1/d)
     963          46 :       k_bac_l            = 0.03_dbl_kind
     964          46 :       algaltype_diatoms  = c0            ! ------------------
     965          46 :       algaltype_sp       = p5            !
     966          46 :       algaltype_phaeo    = p5            !
     967          46 :       nitratetype        = -c1           ! mobility type between
     968          46 :       ammoniumtype       = c1            ! stationary <-->  mobile
     969          46 :       silicatetype       = -c1           !
     970          46 :       dmspptype          = p5            !
     971          46 :       dmspdtype          = -c1           !
     972          46 :       humtype            = c1            !
     973          46 :       doctype_s          = p5            !
     974          46 :       doctype_l          = p5            !
     975          46 :       dontype_protein    = p5            !
     976          46 :       fedtype_1          = p5            !
     977          46 :       feptype_1          = p5            !
     978          46 :       zaerotype_bc1      = c1            !
     979          46 :       zaerotype_bc2      = c1            !
     980          46 :       zaerotype_dust1    = c1            !
     981          46 :       zaerotype_dust2    = c1            !
     982          46 :       zaerotype_dust3    = c1            !
     983          46 :       zaerotype_dust4    = c1            !--------------------
     984          46 :       ratio_C2N_diatoms  = 7.0_dbl_kind  ! algal C to N ratio (mol/mol)
     985          46 :       ratio_C2N_sp       = 7.0_dbl_kind
     986          46 :       ratio_C2N_phaeo    = 7.0_dbl_kind
     987          46 :       ratio_chl2N_diatoms= 2.1_dbl_kind  ! algal chlorophyll to N ratio (mg/mmol)
     988          46 :       ratio_chl2N_sp     = 1.1_dbl_kind
     989          46 :       ratio_chl2N_phaeo  = 0.84_dbl_kind
     990          46 :       F_abs_chl_diatoms  = 2.0_dbl_kind  ! scales absorbed radiation for dEdd
     991          46 :       F_abs_chl_sp       = 4.0_dbl_kind
     992          46 :       F_abs_chl_phaeo    = 5.0_dbl_kind
     993          46 :       ratio_C2N_proteins = 7.0_dbl_kind  ! ratio of C to N in proteins (mol/mol)       
     994             : 
     995             :       !-----------------------------------------------------------------
     996             :       ! read from input file
     997             :       !-----------------------------------------------------------------
     998             : 
     999          46 :       open (nu_nml, file=trim(nml_filename), status='old',iostat=nml_error)
    1000          46 :       if (nml_error /= 0) then
    1001           0 :          nml_error = -1
    1002             :       else
    1003          46 :          nml_error =  1
    1004             :       endif 
    1005             : 
    1006          46 :       print*,'Reading zbgc_nml'
    1007          92 :       do while (nml_error > 0)
    1008          92 :          read(nu_nml, nml=zbgc_nml,iostat=nml_error)
    1009             :       end do
    1010          46 :       if (nml_error == 0) close(nu_nml)
    1011          46 :       if (nml_error /= 0) then
    1012           0 :          print*,'error reading zbgc namelist'
    1013           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1014             :       endif
    1015             : 
    1016             :       !-----------------------------------------------------------------
    1017             :       ! resolve conflicts
    1018             :       !-----------------------------------------------------------------
    1019             :       ! zsalinity and brine
    1020             :       !-----------------------------------------------------------------
    1021          46 :       if (solve_zsal .and. TRZS == 0) then
    1022           0 :          write(nu_diag,*) 'WARNING: solve_zsal=T but 0 zsalinity tracers'
    1023           0 :          write(nu_diag,*) 'WARNING: setting solve_zsal = F'
    1024           0 :          solve_zsal = .false.      
    1025             :       elseif (solve_zsal .and. nblyr < 1)  then
    1026             :          write(nu_diag,*) 'WARNING: solve_zsal=T but 0 zsalinity tracers'
    1027             :          write(nu_diag,*) 'WARNING: setting solve_zsal = F'
    1028             :          solve_zsal = .false.     
    1029             :       endif 
    1030             : 
    1031          46 :       if (solve_zsal .and. ((.not. tr_brine) .or. (ktherm /= 1))) then
    1032           0 :          write(nu_diag,*) 'WARNING: solve_zsal needs tr_brine=T and ktherm=1'
    1033           0 :          write(nu_diag,*) 'WARNING: setting tr_brine=T and ktherm=1'
    1034           0 :          tr_brine = .true.
    1035           0 :          ktherm = 1
    1036             :       endif
    1037             : 
    1038          41 :       if (tr_brine .and. TRBRI == 0 ) then
    1039             :          write(nu_diag,*) &
    1040           0 :             'WARNING: tr_brine=T but no brine height compiled'
    1041             :          write(nu_diag,*) &
    1042           0 :             'WARNING: setting solve_zsal and tr_brine = F'
    1043           0 :          solve_zsal = .false.
    1044           0 :          tr_brine  = .false.
    1045             :       elseif (tr_brine .and. nblyr < 1 ) then
    1046             :          write(nu_diag,*) &
    1047             :             'WARNING: tr_brine=T but no biology layers compiled'
    1048             :          write(nu_diag,*) &
    1049             :             'WARNING: setting solve_zsal and tr_brine = F'
    1050             :          solve_zsal = .false.
    1051             :          tr_brine  = .false.
    1052             :       endif 
    1053             : 
    1054          46 :          write(nu_diag,1010) ' tr_brine                  = ', tr_brine
    1055          46 :       if (tr_brine) then
    1056           5 :          write(nu_diag,1005) ' phi_snow                  = ', phi_snow
    1057             :       endif
    1058          46 :       if (solve_zsal) then
    1059           0 :          write(nu_diag,1010) ' solve_zsal                = ', solve_zsal
    1060           0 :          write(nu_diag,1000) ' grid_oS                   = ', grid_oS
    1061           0 :          write(nu_diag,1005) ' l_skS                     = ', l_skS
    1062             :       endif
    1063             : 
    1064             :       !-----------------------------------------------------------------
    1065             :       ! biogeochemistry
    1066             :       !-----------------------------------------------------------------
    1067             : 
    1068          46 :       if (.not. tr_brine) then
    1069          41 :          if (solve_zbgc) then
    1070           0 :             write(nu_diag,*) 'WARNING: tr_brine = F and solve_zbgc = T'
    1071           0 :             write(nu_diag,*) 'WARNING: setting solve_zbgc = F'
    1072           0 :             solve_zbgc = .false.
    1073             :          endif
    1074          41 :          if (skl_bgc) then
    1075           0 :             write(nu_diag,*) 'WARNING: tr_brine = F and skl_bgc = T'
    1076           0 :             write(nu_diag,*) 'WARNING: setting skl_bgc = F'
    1077           0 :             skl_bgc = .false.
    1078             :          endif
    1079          41 :          if (tr_zaero) then
    1080           0 :             write(nu_diag,*) 'WARNING: tr_brine = F and tr_zaero = T'
    1081           0 :             write(nu_diag,*) 'WARNING: setting tr_zaero = F'
    1082           0 :             tr_zaero = .false.
    1083             :          endif
    1084             :       endif
    1085             : 
    1086          46 :       if ((skl_bgc .AND. solve_zbgc) .or. (skl_bgc .AND. z_tracers)) then
    1087           0 :          print*, 'ERROR: skl_bgc and (solve_zbgc or z_tracers) are both true'
    1088           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1089             :       endif
    1090             : 
    1091          46 :       if (skl_bgc .AND. tr_zaero) then
    1092           0 :          write(nu_diag,*) 'WARNING: skl bgc does not use vertical tracers'
    1093           0 :          write(nu_diag,*) 'WARNING: setting tr_zaero = F'
    1094           0 :          tr_zaero = .false.
    1095             :       endif
    1096             : 
    1097          46 :       if (dEdd_algae .AND. trim(shortwave) /= 'dEdd') then 
    1098           0 :          write(nu_diag,*) 'WARNING: dEdd_algae = T but shortwave /= dEdd'
    1099           0 :          write(nu_diag,*) 'WARNING: setting dEdd_algae = F'
    1100           0 :          dEdd_algae = .false.
    1101             :       endif
    1102             : 
    1103          46 :       if (dEdd_algae .AND. (.NOT. tr_bgc_N) .AND. (.NOT. tr_zaero)) then 
    1104           0 :          write(nu_diag,*) 'WARNING: need tr_bgc_N or tr_zaero for dEdd_algae'
    1105           0 :          write(nu_diag,*) 'WARNING: setting dEdd_algae = F'
    1106           0 :          dEdd_algae = .false.
    1107             :       endif
    1108             : 
    1109          46 :       if (modal_aero .AND. (.NOT. tr_zaero) .AND. (.NOT. tr_aero)) then
    1110           0 :          modal_aero = .false.
    1111             :       endif
    1112             :          
    1113          46 :       if (modal_aero .AND. trim(shortwave) /= 'dEdd') then 
    1114           0 :          write(nu_diag,*) 'WARNING: modal_aero = T but shortwave /= dEdd'
    1115           0 :          write(nu_diag,*) 'WARNING: setting modal_aero = F'
    1116           0 :          modal_aero = .false.
    1117             :       endif
    1118             :       if (n_algae > icepack_max_algae) then
    1119             :          print*, 'error:number of algal types exceeds icepack_max_algae'
    1120             :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1121             :       endif
    1122             :       if (n_doc > icepack_max_doc) then
    1123             :          print*, 'error:number of algal types exceeds icepack_max_doc'
    1124             :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1125             :       endif
    1126             :       if (n_dic > icepack_max_dic) then
    1127             :          print*, 'error:number of dic types exceeds icepack_max_dic'
    1128             :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1129             :       endif
    1130             :       if (n_don > icepack_max_don) then
    1131             :          print*, 'error:number of don types exceeds icepack_max_don'
    1132             :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1133             :       endif
    1134             :       if (n_fed > icepack_max_fe) then
    1135             :          print*, 'error:number of dissolved fe types exceeds icepack_max_fe'
    1136             :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1137             :       endif
    1138             :       if (n_fep > icepack_max_fe) then
    1139             :          print*, 'error:number of particulate fe types exceeds icepack_max_fe'
    1140             :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1141             :       endif
    1142             : 
    1143          45 :       if ((TRBGCS == 0 .and. skl_bgc) .or. (TRALG == 0 .and. skl_bgc)) then
    1144             :          write(nu_diag,*) &
    1145           0 :             'WARNING: skl_bgc=T but 0 bgc or algal tracers compiled'
    1146             :          write(nu_diag,*) &
    1147           0 :             'WARNING: setting skl_bgc = F'
    1148           0 :          skl_bgc = .false.
    1149             :       endif
    1150             : 
    1151          42 :       if ((TRBGCZ == 0 .and. solve_zbgc) .or. (TRALG == 0 .and. solve_zbgc)) then
    1152             :          write(nu_diag,*) &
    1153           0 :             'WARNING: solve_zbgc=T but 0 zbgc or algal tracers compiled'
    1154             :          write(nu_diag,*) &
    1155           0 :             'WARNING: setting solve_zbgc = F'
    1156           0 :          solve_zbgc = .false.
    1157             :       endif
    1158             : 
    1159          46 :       if (solve_zbgc .and. .not. z_tracers) z_tracers = .true.
    1160          46 :       if (skl_bgc .or. solve_zbgc) then
    1161           5 :          tr_bgc_N         = .true.   ! minimum NP biogeochemistry
    1162           5 :          tr_bgc_Nit       = .true.
    1163             :       else
    1164          41 :          tr_bgc_N         = .false.
    1165          41 :          tr_bgc_C         = .false.
    1166          41 :          tr_bgc_chl       = .false.
    1167          41 :          tr_bgc_Nit       = .false.
    1168          41 :          tr_bgc_Am        = .false.
    1169          41 :          tr_bgc_Sil       = .false.
    1170          41 :          tr_bgc_hum       = .false.
    1171          41 :          tr_bgc_DMS       = .false.
    1172          41 :          tr_bgc_PON       = .false.
    1173          41 :          tr_bgc_DON       = .false.
    1174          41 :          tr_bgc_Fe        = .false.
    1175             :       endif
    1176             : 
    1177             :       !-----------------------------------------------------------------
    1178             :       ! z layer aerosols
    1179             :       !-----------------------------------------------------------------
    1180          46 :       if (tr_zaero .and. .not. z_tracers) z_tracers = .true.
    1181             : 
    1182             :       if (n_zaero > icepack_max_aero) then
    1183             :          print*, 'error:number of z aerosols exceeds icepack_max_aero'
    1184             :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1185             :       endif         
    1186             : 
    1187             :       if (skl_bgc .and. n_bgc < 2) then
    1188             :          write (nu_diag,*) ' '
    1189             :          write (nu_diag,*) 'comp_ice must have number of bgc tracers >= 2'
    1190             :          write (nu_diag,*) 'number of bgc tracers compiled:',n_bgc
    1191             :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1192             :       endif
    1193             : 
    1194             :       if (solve_zbgc .and. n_bgc < 2) then
    1195             :          write (nu_diag,*) ' '
    1196             :          write (nu_diag,*) 'comp_ice must have number of zbgc tracers >= 2'
    1197             :          write (nu_diag,*) 'number of bgc tracers compiled:',n_bgc
    1198             :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1199             :       endif
    1200             : 
    1201          43 :       if (tr_zaero .and. TRZAERO <  1) then
    1202           0 :          write (nu_diag,*) ' '
    1203           0 :          write (nu_diag,*) 'comp_ice must have number of TRZAERO > 0'
    1204           0 :          write (nu_diag,*) 'in order to solve z aerosols:',TRZAERO
    1205           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1206             :       endif
    1207             : 
    1208             :       !-----------------------------------------------------------------
    1209             :       ! end conflict resolution
    1210             :       !-----------------------------------------------------------------
    1211             :       !-----------------------------------------------------------------
    1212             :       ! set Icepack values
    1213             :       !-----------------------------------------------------------------
    1214             : 
    1215             :       call icepack_init_parameters(ktherm_in=ktherm, shortwave_in=shortwave, &
    1216             :            scale_bgc_in=scale_bgc, skl_bgc_in=skl_bgc, z_tracers_in=z_tracers, &
    1217             :            dEdd_algae_in=dEdd_algae, solve_zbgc_in=solve_zbgc, &
    1218             :            bgc_flux_type_in=bgc_flux_type, grid_o_in=grid_o, l_sk_in=l_sk, &
    1219             :            initbio_frac_in=initbio_frac, frazil_scav_in=frazil_scav, &
    1220             :            grid_oS_in=grid_oS, l_skS_in=l_skS, phi_snow_in=phi_snow, &
    1221             :            algal_vel_in=algal_vel, R_dFe2dust_in=R_dFe2dust, &
    1222             :            dustFe_sol_in=dustFe_sol, T_max_in=T_max, fsal_in=fsal, &
    1223             :            op_dep_min_in=op_dep_min, fr_graze_s_in=fr_graze_s, &
    1224             :            fr_graze_e_in=fr_graze_e, fr_mort2min_in=fr_mort2min, &
    1225             :            fr_dFe_in=fr_dFe, k_nitrif_in=k_nitrif, t_iron_conv_in=t_iron_conv, &
    1226             :            max_loss_in=max_loss, max_dfe_doc1_in=max_dfe_doc1, fr_resp_in=fr_resp, &
    1227             :            fr_resp_s_in=fr_resp_s, y_sk_DMS_in=y_sk_DMS, t_sk_conv_in=t_sk_conv, &
    1228          46 :            t_sk_ox_in=t_sk_ox, modal_aero_in=modal_aero)
    1229          46 :       call icepack_warnings_flush(nu_diag)
    1230          46 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
    1231           0 :           file=__FILE__, line=__LINE__)
    1232             : 
    1233             :       !-----------------------------------------------------------------
    1234             :       ! initialize zbgc tracer indices
    1235             :       !----------------------------------------------------------------- 
    1236             : 
    1237          46 :       ntrcr_o = ntrcr
    1238          46 :       nt_fbri = 0
    1239          46 :       if (tr_brine) then
    1240           5 :           nt_fbri = ntrcr + 1   ! ice volume fraction with salt
    1241           5 :           ntrcr = ntrcr + 1
    1242           5 :           trcr_depend(nt_fbri)   = 1   ! volume-weighted
    1243           5 :           trcr_base  (nt_fbri,1) = c0  ! volume-weighted
    1244           5 :           trcr_base  (nt_fbri,2) = c1  ! volume-weighted
    1245           5 :           trcr_base  (nt_fbri,3) = c0  ! volume-weighted
    1246           5 :           n_trcr_strata(nt_fbri) = 0
    1247           5 :           nt_strata  (nt_fbri,1) = 0
    1248           5 :           nt_strata  (nt_fbri,2) = 0
    1249             :       endif
    1250             : 
    1251          46 :       ntd = 0                    ! if nt_fbri /= 0 then use fbri dependency
    1252          46 :       if (nt_fbri == 0) ntd = -1 ! otherwise make tracers depend on ice volume
    1253             : 
    1254          46 :       nt_bgc_S = 0
    1255          46 :       if (solve_zsal) then       ! .true. only if tr_brine = .true.
    1256           0 :           nt_bgc_S = ntrcr + 1
    1257           0 :           ntrcr = ntrcr + nblyr
    1258           0 :           do k = 1,nblyr
    1259           0 :              trcr_depend(nt_bgc_S + k - 1) = 2 + nt_fbri + ntd
    1260           0 :              trcr_base  (nt_bgc_S,1) = c0  ! default: ice area
    1261           0 :              trcr_base  (nt_bgc_S,2) = c1 
    1262           0 :              trcr_base  (nt_bgc_S,3) = c0  
    1263           0 :              n_trcr_strata(nt_bgc_S) = 1
    1264           0 :              nt_strata(nt_bgc_S,1) = nt_fbri
    1265           0 :              nt_strata(nt_bgc_S,2) = 0
    1266             :           enddo
    1267             :       endif 
    1268             : 
    1269             :       !-----------------------------------------------------------------
    1270             :       ! biogeochemistry
    1271             :       !-----------------------------------------------------------------
    1272             : 
    1273          46 :       nbtrcr = 0
    1274          46 :       nbtrcr_sw = 0
    1275             : 
    1276             :       ! vectors of size icepack_max_algae
    1277          46 :       nlt_bgc_N(:) = 0
    1278          46 :       nlt_bgc_C(:) = 0
    1279          46 :       nlt_bgc_chl(:) = 0
    1280          46 :       nt_bgc_N(:) = 0
    1281          46 :       nt_bgc_C(:) = 0
    1282          46 :       nt_bgc_chl(:) = 0
    1283             : 
    1284             :       ! vectors of size icepack_max_dic
    1285          46 :       nlt_bgc_DIC(:) = 0
    1286          46 :       nt_bgc_DIC(:) = 0
    1287             : 
    1288             :       ! vectors of size icepack_max_doc
    1289          46 :       nlt_bgc_DOC(:) = 0
    1290          46 :       nt_bgc_DOC(:) = 0
    1291             : 
    1292             :       ! vectors of size icepack_max_don
    1293          46 :       nlt_bgc_DON(:) = 0
    1294          46 :       nt_bgc_DON(:) = 0
    1295             : 
    1296             :       ! vectors of size icepack_max_fe 
    1297          46 :       nlt_bgc_Fed(:) = 0
    1298          46 :       nlt_bgc_Fep(:) = 0
    1299          46 :       nt_bgc_Fed(:) = 0
    1300          46 :       nt_bgc_Fep(:) = 0
    1301             : 
    1302             :       ! vectors of size icepack_max_aero
    1303          46 :       nlt_zaero(:) = 0
    1304          46 :       nlt_zaero_sw(:) = 0
    1305          46 :       nt_zaero(:) = 0
    1306          46 :       nt_zaero_sw(:) = 0
    1307             : 
    1308          46 :       nlt_bgc_Nit    = 0
    1309          46 :       nlt_bgc_Am     = 0
    1310          46 :       nlt_bgc_Sil    = 0
    1311          46 :       nlt_bgc_DMSPp  = 0
    1312          46 :       nlt_bgc_DMSPd  = 0
    1313          46 :       nlt_bgc_DMS    = 0
    1314          46 :       nlt_bgc_PON    = 0
    1315          46 :       nlt_bgc_hum    = 0
    1316          46 :       nlt_chl_sw     = 0
    1317          46 :       bio_index(:)   = 0
    1318          46 :       bio_index_o(:) = 0
    1319             : 
    1320          46 :       nt_bgc_Nit    = 0
    1321          46 :       nt_bgc_Am     = 0
    1322          46 :       nt_bgc_Sil    = 0
    1323          46 :       nt_bgc_DMSPp  = 0
    1324          46 :       nt_bgc_DMSPd  = 0
    1325          46 :       nt_bgc_DMS    = 0
    1326          46 :       nt_bgc_PON    = 0
    1327          46 :       nt_bgc_hum    = 0
    1328          46 :       nt_zbgc_frac  = 0
    1329             : 
    1330             :       !-----------------------------------------------------------------
    1331             :       ! Define array parameters
    1332             :       !-----------------------------------------------------------------
    1333             : 
    1334          46 :       R_Si2N(1) = ratio_Si2N_diatoms
    1335          46 :       R_Si2N(2) = ratio_Si2N_sp
    1336          46 :       R_Si2N(3) = ratio_Si2N_phaeo
    1337             : 
    1338          46 :       R_S2N(1) = ratio_S2N_diatoms
    1339          46 :       R_S2N(2) = ratio_S2N_sp
    1340          46 :       R_S2N(3) = ratio_S2N_phaeo
    1341             : 
    1342          46 :       R_Fe2C(1) = ratio_Fe2C_diatoms
    1343          46 :       R_Fe2C(2) = ratio_Fe2C_sp
    1344          46 :       R_Fe2C(3) = ratio_Fe2C_phaeo
    1345             : 
    1346          46 :       R_Fe2N(1) = ratio_Fe2N_diatoms
    1347          46 :       R_Fe2N(2) = ratio_Fe2N_sp
    1348          46 :       R_Fe2N(3) = ratio_Fe2N_phaeo
    1349             : 
    1350          46 :       R_C2N(1) = ratio_C2N_diatoms
    1351          46 :       R_C2N(2) = ratio_C2N_sp
    1352          46 :       R_C2N(3) = ratio_C2N_phaeo
    1353             : 
    1354          46 :       R_chl2N(1) = ratio_chl2N_diatoms
    1355          46 :       R_chl2N(2) = ratio_chl2N_sp
    1356          46 :       R_chl2N(3) = ratio_chl2N_phaeo
    1357             : 
    1358          46 :       F_abs_chl(1) = F_abs_chl_diatoms
    1359          46 :       F_abs_chl(2) = F_abs_chl_sp
    1360          46 :       F_abs_chl(3) = F_abs_chl_phaeo
    1361             : 
    1362          46 :       R_Fe2DON(1) = ratio_Fe2DON
    1363          46 :       R_C2N_DON(1) = ratio_C2N_proteins
    1364             :      
    1365          46 :       R_Fe2DOC(1) = ratio_Fe2DOC_s
    1366          46 :       R_Fe2DOC(2) = ratio_Fe2DOC_l
    1367          46 :       R_Fe2DOC(3) = c0
    1368             : 
    1369          46 :       chlabs(1) = chlabs_diatoms
    1370          46 :       chlabs(2) = chlabs_sp
    1371          46 :       chlabs(3) = chlabs_phaeo
    1372             : 
    1373          46 :       alpha2max_low(1) = alpha2max_low_diatoms
    1374          46 :       alpha2max_low(2) = alpha2max_low_sp
    1375          46 :       alpha2max_low(3) = alpha2max_low_phaeo
    1376             : 
    1377          46 :       beta2max(1) = beta2max_diatoms
    1378          46 :       beta2max(2) = beta2max_sp
    1379          46 :       beta2max(3) = beta2max_phaeo
    1380             : 
    1381          46 :       mu_max(1) = mu_max_diatoms
    1382          46 :       mu_max(2) = mu_max_sp
    1383          46 :       mu_max(3) = mu_max_phaeo
    1384             : 
    1385          46 :       grow_Tdep(1) = grow_Tdep_diatoms
    1386          46 :       grow_Tdep(2) = grow_Tdep_sp
    1387          46 :       grow_Tdep(3) = grow_Tdep_phaeo
    1388             : 
    1389          46 :       fr_graze(1) = fr_graze_diatoms
    1390          46 :       fr_graze(2) = fr_graze_sp
    1391          46 :       fr_graze(3) = fr_graze_phaeo
    1392             : 
    1393          46 :       mort_pre(1) = mort_pre_diatoms
    1394          46 :       mort_pre(2) = mort_pre_sp
    1395          46 :       mort_pre(3) = mort_pre_phaeo
    1396             : 
    1397          46 :       mort_Tdep(1) = mort_Tdep_diatoms
    1398          46 :       mort_Tdep(2) = mort_Tdep_sp
    1399          46 :       mort_Tdep(3) = mort_Tdep_phaeo
    1400             : 
    1401          46 :       k_exude(1) = k_exude_diatoms
    1402          46 :       k_exude(2) = k_exude_sp
    1403          46 :       k_exude(3) = k_exude_phaeo
    1404             : 
    1405          46 :       K_Nit(1) = K_Nit_diatoms
    1406          46 :       K_Nit(2) = K_Nit_sp
    1407          46 :       K_Nit(3) = K_Nit_phaeo
    1408             : 
    1409          46 :       K_Am(1) = K_Am_diatoms
    1410          46 :       K_Am(2) = K_Am_sp
    1411          46 :       K_Am(3) = K_Am_phaeo
    1412             : 
    1413          46 :       K_Sil(1) = K_Sil_diatoms
    1414          46 :       K_Sil(2) = K_Sil_sp
    1415          46 :       K_Sil(3) = K_Sil_phaeo
    1416             : 
    1417          46 :       K_Fe(1) = K_Fe_diatoms
    1418          46 :       K_Fe(2) = K_Fe_sp
    1419          46 :       K_Fe(3) = K_Fe_phaeo
    1420             : 
    1421          46 :       f_don(1) = f_don_protein
    1422          46 :       kn_bac(1) = kn_bac_protein
    1423          46 :       f_don_Am(1) = f_don_Am_protein
    1424             : 
    1425          46 :       f_exude(1) = f_exude_s
    1426          46 :       f_exude(2) = f_exude_l
    1427          46 :       k_bac(1) = k_bac_s
    1428          46 :       k_bac(2) = k_bac_l
    1429             : 
    1430          92 :       dictype(:) = -c1
    1431             :       
    1432          46 :       algaltype(1) = algaltype_diatoms
    1433          46 :       algaltype(2) = algaltype_sp
    1434          46 :       algaltype(3) = algaltype_phaeo
    1435             : 
    1436          46 :       doctype(1) = doctype_s
    1437          46 :       doctype(2) = doctype_l
    1438             :  
    1439          46 :       dontype(1) = dontype_protein
    1440             : 
    1441          46 :       fedtype(1) = fedtype_1
    1442          46 :       feptype(1) = feptype_1
    1443             : 
    1444          46 :       zaerotype(1) = zaerotype_bc1
    1445          46 :       zaerotype(2) = zaerotype_bc2
    1446          46 :       zaerotype(3) = zaerotype_dust1
    1447          46 :       zaerotype(4) = zaerotype_dust2
    1448          46 :       zaerotype(5) = zaerotype_dust3
    1449          46 :       zaerotype(6) = zaerotype_dust4
    1450             : 
    1451             :       call icepack_init_zbgc ( &
    1452             :            R_Si2N_in=R_Si2N, &
    1453             :            R_S2N_in=R_S2N, R_Fe2C_in=R_Fe2C, R_Fe2N_in=R_Fe2N, R_C2N_in=R_C2N, &
    1454             :            R_chl2N_in=R_chl2N, F_abs_chl_in=F_abs_chl, R_Fe2DON_in=R_Fe2DON, &
    1455             :            R_C2N_DON_in=R_C2N_DON, &
    1456             :            R_Fe2DOC_in=R_Fe2DOC, &
    1457             :            chlabs_in=chlabs, alpha2max_low_in=alpha2max_low, beta2max_in=beta2max, &
    1458             :            mu_max_in=mu_max, grow_Tdep_in=grow_Tdep, fr_graze_in=fr_graze, &
    1459             :            mort_pre_in=mort_pre, &
    1460             :            mort_Tdep_in=mort_Tdep, k_exude_in=k_exude, &
    1461             :            K_Nit_in=K_Nit, K_Am_in=K_Am, K_sil_in=K_Sil, K_Fe_in=K_Fe, &
    1462             :            f_don_in=f_don, kn_bac_in=kn_bac, f_don_Am_in=f_don_Am, f_exude_in=f_exude, &
    1463             :            k_bac_in=k_bac, &
    1464             :            fr_resp_in=fr_resp, algal_vel_in=algal_vel, R_dFe2dust_in=R_dFe2dust, &
    1465             :            dustFe_sol_in=dustFe_sol, T_max_in=T_max, fr_mort2min_in=fr_mort2min, &
    1466             :            fr_dFe_in=fr_dFe, op_dep_min_in=op_dep_min, &
    1467             :            fr_graze_s_in=fr_graze_s, fr_graze_e_in=fr_graze_e, &
    1468             :            k_nitrif_in=k_nitrif, t_iron_conv_in=t_iron_conv, &
    1469             :            max_loss_in=max_loss, max_dfe_doc1_in=max_dfe_doc1, &
    1470             :            fr_resp_s_in=fr_resp_s, y_sk_DMS_in=y_sk_DMS, &
    1471          46 :            t_sk_conv_in=t_sk_conv, t_sk_ox_in=t_sk_ox, fsal_in=fsal)
    1472          46 :       call icepack_warnings_flush(nu_diag)
    1473          46 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
    1474           0 :           file=__FILE__, line=__LINE__)
    1475             : 
    1476          46 :       if (skl_bgc) then
    1477             : 
    1478           1 :          nk = 1
    1479           1 :          nt_depend = 0
    1480             : 
    1481           1 :          if (dEdd_algae) then 
    1482           0 :            nlt_chl_sw = 1
    1483           0 :            nbtrcr_sw = nilyr+nslyr+2  ! only the bottom layer will be nonzero
    1484             :          endif  
    1485             :          
    1486          45 :       elseif (z_tracers) then ! defined on nblyr+1 in ice 
    1487             :                               ! and 2 snow layers (snow surface + interior)
    1488             : 
    1489           4 :          nk = nblyr + 1
    1490           4 :          nt_depend = 2 + nt_fbri + ntd 
    1491             : 
    1492           4 :          if (tr_bgc_N) then
    1493           4 :             if (dEdd_algae) then
    1494           0 :                nlt_chl_sw = 1
    1495           0 :                nbtrcr_sw =  nilyr+nslyr+2
    1496             :             endif
    1497             :          endif ! tr_bgc_N
    1498             : 
    1499             :       endif ! skl_bgc or z_tracers
    1500             : 
    1501          46 :       if (skl_bgc .or. z_tracers) then
    1502             : 
    1503             :       !-----------------------------------------------------------------
    1504             :       ! assign tracer indices and dependencies
    1505             :       ! bgc_tracer_type: < 0  purely mobile , >= 0 stationary 
    1506             :       !------------------------------------------------------------------
    1507             : 
    1508           5 :       if (tr_bgc_N) then
    1509          20 :          do mm = 1, n_algae      
    1510             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1511           9 :                                nt_bgc_N(mm),    nlt_bgc_N(mm), &
    1512           9 :                                algaltype(mm),   nt_depend,     &
    1513             :                                ntrcr,           nbtrcr,        &
    1514             :                                bgc_tracer_type, trcr_depend,   &
    1515             :                                trcr_base,       n_trcr_strata, &
    1516          15 :                                nt_strata,       bio_index)
    1517          20 :             bio_index_o(nlt_bgc_N(mm)) = mm
    1518             :          enddo   ! mm
    1519             :       endif ! tr_bgc_N
    1520             : 
    1521           5 :       if (tr_bgc_Nit) then
    1522             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1523             :                                nt_bgc_Nit,      nlt_bgc_Nit,   &
    1524             :                                nitratetype,     nt_depend,     &
    1525             :                                ntrcr,           nbtrcr,        &
    1526             :                                bgc_tracer_type, trcr_depend,   &
    1527             :                                trcr_base,       n_trcr_strata, &
    1528           5 :                                nt_strata,       bio_index)
    1529           5 :             bio_index_o(nlt_bgc_Nit) = icepack_max_algae + 1
    1530             :       endif ! tr_bgc_Nit
    1531             :          
    1532           5 :       if (tr_bgc_C) then
    1533             :        !
    1534             :        ! Algal C is not yet distinct from algal N
    1535             :        ! * Reqires exudation and/or changing C:N ratios
    1536             :        ! for implementation
    1537             :        !
    1538             :        !  do mm = 1,n_algae      
    1539             :        !     call init_bgc_trcr(nk,              nt_fbri,       &
    1540             :        !                        nt_bgc_C(mm),    nlt_bgc_C(mm), &
    1541             :        !                        algaltype(mm),   nt_depend,     &
    1542             :        !                        ntrcr,           nbtrcr,        &
    1543             :        !                        bgc_tracer_type, trcr_depend,   &
    1544             :        !                        trcr_base,       n_trcr_strata, &
    1545             :        !                        nt_strata,       bio_index)
    1546             :        !     bio_index_o(nlt_bgc_C(mm)) = icepack_max_algae + 1 + mm
    1547             :        !  enddo   ! mm
    1548             : 
    1549          15 :          do mm = 1, n_doc
    1550             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1551           6 :                                nt_bgc_DOC(mm),  nlt_bgc_DOC(mm), &
    1552           6 :                                doctype(mm),     nt_depend,     &
    1553             :                                ntrcr,           nbtrcr,        &
    1554             :                                bgc_tracer_type, trcr_depend,   &
    1555             :                                trcr_base,       n_trcr_strata, &
    1556          10 :                                nt_strata,       bio_index)
    1557          15 :             bio_index_o(nlt_bgc_DOC(mm)) = icepack_max_algae + 1 + mm
    1558             :          enddo   ! mm
    1559           5 :          do mm = 1, n_dic
    1560             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1561           0 :                                nt_bgc_DIC(mm),  nlt_bgc_DIC(mm), &
    1562           0 :                                dictype(mm),     nt_depend,     &
    1563             :                                ntrcr,           nbtrcr,        &
    1564             :                                bgc_tracer_type, trcr_depend,   &
    1565             :                                trcr_base,       n_trcr_strata, &
    1566           0 :                                nt_strata,       bio_index)
    1567           5 :             bio_index_o(nlt_bgc_DIC(mm)) = icepack_max_algae + icepack_max_doc + 1 + mm
    1568             :          enddo   ! mm
    1569             :       endif      ! tr_bgc_C
    1570             : 
    1571           5 :       if (tr_bgc_chl) then
    1572           0 :          do mm = 1, n_algae
    1573             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1574           0 :                                nt_bgc_chl(mm),  nlt_bgc_chl(mm), &
    1575           0 :                                algaltype(mm),   nt_depend,     &
    1576             :                                ntrcr,           nbtrcr,        &
    1577             :                                bgc_tracer_type, trcr_depend,   &
    1578             :                                trcr_base,       n_trcr_strata, &
    1579           0 :                                nt_strata,       bio_index)
    1580           0 :             bio_index_o(nlt_bgc_chl(mm)) = icepack_max_algae + 1 + icepack_max_doc + icepack_max_dic + mm
    1581             :          enddo   ! mm
    1582             :       endif      ! tr_bgc_chl
    1583             : 
    1584           5 :       if (tr_bgc_Am) then
    1585             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1586             :                                nt_bgc_Am,       nlt_bgc_Am,    &
    1587             :                                ammoniumtype,    nt_depend,     &
    1588             :                                ntrcr,           nbtrcr,        &
    1589             :                                bgc_tracer_type, trcr_depend,   &
    1590             :                                trcr_base,       n_trcr_strata, &
    1591           5 :                                nt_strata,       bio_index)
    1592           5 :             bio_index_o(nlt_bgc_Am) = 2*icepack_max_algae + icepack_max_doc + icepack_max_dic + 2
    1593             :       endif    
    1594           5 :       if (tr_bgc_Sil) then
    1595             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1596             :                                nt_bgc_Sil,      nlt_bgc_Sil,   &
    1597             :                                silicatetype,    nt_depend,     &
    1598             :                                ntrcr,           nbtrcr,        &
    1599             :                                bgc_tracer_type, trcr_depend,   &
    1600             :                                trcr_base,       n_trcr_strata, &
    1601           5 :                                nt_strata,       bio_index)
    1602           5 :             bio_index_o(nlt_bgc_Sil) = 2*icepack_max_algae + icepack_max_doc + icepack_max_dic + 3
    1603             :       endif    
    1604           5 :       if (tr_bgc_DMS) then   ! all together
    1605             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1606             :                                nt_bgc_DMSPp,    nlt_bgc_DMSPp, &
    1607             :                                dmspptype,       nt_depend,     &
    1608             :                                ntrcr,           nbtrcr,        &
    1609             :                                bgc_tracer_type, trcr_depend,   &
    1610             :                                trcr_base,       n_trcr_strata, &
    1611           5 :                                nt_strata,       bio_index)
    1612           5 :             bio_index_o(nlt_bgc_DMSPp) = 2*icepack_max_algae + icepack_max_doc + icepack_max_dic + 4
    1613             : 
    1614             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1615             :                                nt_bgc_DMSPd,    nlt_bgc_DMSPd, &
    1616             :                                dmspdtype,       nt_depend,     &
    1617             :                                ntrcr,           nbtrcr,        &
    1618             :                                bgc_tracer_type, trcr_depend,   &
    1619             :                                trcr_base,       n_trcr_strata, &
    1620           5 :                                nt_strata,       bio_index)
    1621           5 :             bio_index_o(nlt_bgc_DMSPd) = 2*icepack_max_algae + icepack_max_doc + icepack_max_dic + 5
    1622             : 
    1623             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1624             :                                nt_bgc_DMS,      nlt_bgc_DMS,   &
    1625             :                                dmspdtype,       nt_depend,     &
    1626             :                                ntrcr,           nbtrcr,        &
    1627             :                                bgc_tracer_type, trcr_depend,   &
    1628             :                                trcr_base,       n_trcr_strata, &
    1629           5 :                                nt_strata,       bio_index)
    1630           5 :             bio_index_o(nlt_bgc_DMS) = 2*icepack_max_algae + icepack_max_doc + icepack_max_dic + 6
    1631             :       endif    
    1632           5 :       if (tr_bgc_PON) then
    1633             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1634             :                                nt_bgc_PON,      nlt_bgc_PON, &
    1635             :                                nitratetype,     nt_depend,     &
    1636             :                                ntrcr,           nbtrcr,        &
    1637             :                                bgc_tracer_type, trcr_depend,   &
    1638             :                                trcr_base,       n_trcr_strata, &
    1639           5 :                                nt_strata,       bio_index)
    1640           5 :             bio_index_o(nlt_bgc_PON) =  2*icepack_max_algae + icepack_max_doc + icepack_max_dic + 7
    1641             :       endif
    1642           5 :       if (tr_bgc_DON) then
    1643          10 :          do mm = 1, n_don
    1644             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1645           3 :                                nt_bgc_DON(mm),  nlt_bgc_DON(mm), &
    1646           3 :                                dontype(mm),     nt_depend,     &
    1647             :                                ntrcr,           nbtrcr,        &
    1648             :                                bgc_tracer_type, trcr_depend,   &
    1649             :                                trcr_base,       n_trcr_strata, &
    1650           5 :                                nt_strata,       bio_index)
    1651          10 :             bio_index_o(nlt_bgc_DON(mm)) = 2*icepack_max_algae + icepack_max_doc + icepack_max_dic + 7 + mm
    1652             :          enddo   ! mm
    1653             :       endif      ! tr_bgc_DON
    1654           5 :       if (tr_bgc_Fe) then
    1655           4 :          do mm = 1, n_fed
    1656             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1657           2 :                                nt_bgc_Fed(mm),  nlt_bgc_Fed(mm), &
    1658           2 :                                fedtype(mm),     nt_depend,     &
    1659             :                                ntrcr,           nbtrcr,        &
    1660             :                                bgc_tracer_type, trcr_depend,   &
    1661             :                                trcr_base,       n_trcr_strata, &
    1662           2 :                                nt_strata,       bio_index)
    1663           4 :             bio_index_o(nlt_bgc_Fed(mm)) = 2*icepack_max_algae + icepack_max_doc + icepack_max_dic &
    1664           6 :                                          + icepack_max_don + 7 + mm
    1665             :          enddo   ! mm
    1666           4 :          do mm = 1, n_fep
    1667             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1668           2 :                                nt_bgc_Fep(mm),  nlt_bgc_Fep(mm), &
    1669           2 :                                feptype(mm),     nt_depend,     &
    1670             :                                ntrcr,           nbtrcr,        &
    1671             :                                bgc_tracer_type, trcr_depend,   &
    1672             :                                trcr_base,       n_trcr_strata, &
    1673           2 :                                nt_strata,       bio_index)
    1674           4 :             bio_index_o(nlt_bgc_Fep(mm)) = 2*icepack_max_algae + icepack_max_doc + icepack_max_dic &
    1675           6 :                                          + icepack_max_don + icepack_max_fe + 7 + mm
    1676             :          enddo   ! mm
    1677             :       endif      ! tr_bgc_Fe 
    1678             :   
    1679           5 :       if (tr_bgc_hum) then
    1680             :             call init_bgc_trcr(nk,              nt_fbri,       &
    1681             :                                nt_bgc_hum,      nlt_bgc_hum,   &
    1682             :                                humtype,         nt_depend,     &
    1683             :                                ntrcr,           nbtrcr,        &
    1684             :                                bgc_tracer_type, trcr_depend,   &
    1685             :                                trcr_base,       n_trcr_strata, &
    1686           5 :                                nt_strata,       bio_index)
    1687           3 :             bio_index_o(nlt_bgc_hum) =   2*icepack_max_algae + icepack_max_doc + 8 + icepack_max_dic &
    1688           5 :                                          + icepack_max_don + 2*icepack_max_fe + icepack_max_aero 
    1689             :       endif
    1690             :       endif  ! skl_bgc or z_tracers
    1691             : 
    1692          46 :       if (z_tracers) then ! defined on nblyr+1 in ice 
    1693             :                               ! and 2 snow layers (snow surface + interior)
    1694             : 
    1695           4 :          nk = nblyr + 1
    1696           4 :          nt_depend = 2 + nt_fbri + ntd 
    1697             : 
    1698             :          ! z layer aerosols
    1699           4 :          if (tr_zaero) then
    1700          21 :             do mm = 1, n_zaero
    1701          18 :                if (dEdd_algae) then
    1702           0 :                   nlt_zaero_sw(mm) = nbtrcr_sw + 1
    1703           0 :                   nbtrcr_sw = nbtrcr_sw + nilyr + nslyr+2
    1704             :                endif
    1705             :                call init_bgc_trcr(nk,              nt_fbri,       &
    1706           6 :                                   nt_zaero(mm),    nlt_zaero(mm), &
    1707           6 :                                   zaerotype(mm),   nt_depend,     &
    1708             :                                   ntrcr,           nbtrcr,        &
    1709             :                                   bgc_tracer_type, trcr_depend,   &
    1710             :                                   trcr_base,       n_trcr_strata, &
    1711          18 :                                   nt_strata,       bio_index)
    1712          12 :                bio_index_o(nlt_zaero(mm)) = 2*icepack_max_algae + icepack_max_doc + icepack_max_dic &
    1713          27 :                                           + icepack_max_don + 2*icepack_max_fe + 7 + mm
    1714             :             enddo   ! mm
    1715             :          endif      ! tr_zaero
    1716             : 
    1717             : !echmod keep trcr indices etc here but move zbgc_frac_init, zbgc_init_frac, tau_ret, tau_rel to icepack
    1718           4 :          if (nbtrcr > 0) then
    1719           4 :             nt_zbgc_frac = ntrcr + 1
    1720           4 :             ntrcr = ntrcr + nbtrcr
    1721          80 :             do k = 1,nbtrcr 
    1722          76 :                zbgc_frac_init(k) = c1   
    1723          76 :                trcr_depend(nt_zbgc_frac+k-1) =  2+nt_fbri 
    1724          76 :                trcr_base(nt_zbgc_frac+ k - 1,1)  = c0
    1725          76 :                trcr_base(nt_zbgc_frac+ k - 1,2)  = c1
    1726          76 :                trcr_base(nt_zbgc_frac+ k - 1,3)  = c0
    1727          76 :                n_trcr_strata(nt_zbgc_frac+ k - 1)= 1  
    1728          76 :                nt_strata(nt_zbgc_frac+ k - 1,1)  = nt_fbri
    1729          76 :                nt_strata(nt_zbgc_frac+ k - 1,2)  = 0 
    1730          76 :                tau_ret(k) = c1
    1731          76 :                tau_rel(k) = c1
    1732          80 :                if (bgc_tracer_type(k) >=  c0 .and. bgc_tracer_type(k) < p5) then
    1733           4 :                   tau_ret(k) = tau_min
    1734           4 :                   tau_rel(k) = tau_max
    1735           4 :                   zbgc_frac_init(k) = c1
    1736          72 :                elseif (bgc_tracer_type(k) >= p5 .and. bgc_tracer_type(k) < c1) then
    1737          26 :                   tau_ret(k) = tau_min
    1738          26 :                   tau_rel(k) = tau_min
    1739          26 :                   zbgc_frac_init(k) = c1
    1740          46 :                elseif (bgc_tracer_type(k) >= c1 .and. bgc_tracer_type(k) < c2) then
    1741          26 :                   tau_ret(k) = tau_max
    1742          26 :                   tau_rel(k) = tau_min
    1743          26 :                   zbgc_frac_init(k) = c1
    1744          20 :                elseif (bgc_tracer_type(k) >= c2 ) then
    1745           0 :                   tau_ret(k) = tau_max
    1746           0 :                   tau_rel(k) = tau_max
    1747           0 :                   zbgc_frac_init(k) = c1
    1748             :                endif
    1749             :             enddo
    1750             :          endif
    1751             : 
    1752             :       endif ! z_tracers
    1753             : 
    1754         138 :       do k = 1, nbtrcr
    1755          92 :          zbgc_init_frac(k) = frazil_scav
    1756         138 :          if (bgc_tracer_type(k) < c0)  zbgc_init_frac(k) = initbio_frac
    1757             :       enddo  
    1758             : 
    1759             :       !-----------------------------------------------------------------
    1760             :       ! set values in icepack
    1761             :       !-----------------------------------------------------------------
    1762             : 
    1763             :       call icepack_init_zbgc( &
    1764             :            zbgc_init_frac_in=zbgc_init_frac, tau_ret_in=tau_ret, tau_rel_in=tau_rel, &
    1765          46 :            zbgc_frac_init_in=zbgc_frac_init, bgc_tracer_type_in=bgc_tracer_type)
    1766          46 :       call icepack_warnings_flush(nu_diag)
    1767          46 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
    1768           0 :           file=__FILE__, line=__LINE__)
    1769             : 
    1770             :       call icepack_init_tracer_sizes( &
    1771             :            n_algae_in=n_algae,                                                                   &
    1772             :            n_DOC_in=n_DOC,             n_DON_in=n_DON,               n_DIC_in=n_DIC,             &
    1773             :            n_fed_in=n_fed,             n_fep_in=n_fep,               n_zaero_in=n_zaero,         &
    1774          46 :            ntrcr_in=ntrcr, ntrcr_o_in=ntrcr_o, nbtrcr_in=nbtrcr, nbtrcr_sw_in=nbtrcr_sw)
    1775          46 :       call icepack_warnings_flush(nu_diag)
    1776          46 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
    1777           0 :           file=__FILE__, line=__LINE__)
    1778             : 
    1779             :       call icepack_init_tracer_flags( &
    1780             :            tr_brine_in  =tr_brine, &
    1781             :            tr_bgc_Nit_in=tr_bgc_Nit, tr_bgc_Am_in =tr_bgc_Am,  tr_bgc_Sil_in=tr_bgc_Sil,   &
    1782             :            tr_bgc_DMS_in=tr_bgc_DMS, tr_bgc_PON_in=tr_bgc_PON,                             &
    1783             :            tr_bgc_N_in  =tr_bgc_N,   tr_bgc_C_in  =tr_bgc_C,   tr_bgc_chl_in=tr_bgc_chl,   &
    1784             :            tr_bgc_DON_in=tr_bgc_DON, tr_bgc_Fe_in =tr_bgc_Fe,  tr_zaero_in  =tr_zaero,     &
    1785          46 :            tr_bgc_hum_in=tr_bgc_hum)
    1786          46 :       call icepack_warnings_flush(nu_diag)
    1787          46 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
    1788           0 :           file=__FILE__, line=__LINE__)
    1789             : 
    1790             :       call icepack_init_tracer_indices( &
    1791             :            nt_fbri_in=nt_fbri,                                                                   &  
    1792             :            nt_bgc_Nit_in=nt_bgc_Nit,   nt_bgc_Am_in=nt_bgc_Am,       nt_bgc_Sil_in=nt_bgc_Sil,   &
    1793             :            nt_bgc_DMS_in=nt_bgc_DMS,   nt_bgc_PON_in=nt_bgc_PON,     nt_bgc_S_in=nt_bgc_S,       &
    1794             :            nt_bgc_N_in=nt_bgc_N,       nt_bgc_chl_in=nt_bgc_chl,     nt_bgc_hum_in=nt_bgc_hum,   &
    1795             :            nt_bgc_DOC_in=nt_bgc_DOC,   nt_bgc_DON_in=nt_bgc_DON,     nt_bgc_DIC_in=nt_bgc_DIC,   &
    1796             :            nt_zaero_in=nt_zaero,       nt_bgc_DMSPp_in=nt_bgc_DMSPp, nt_bgc_DMSPd_in=nt_bgc_DMSPd, &
    1797             :            nt_bgc_Fed_in=nt_bgc_Fed,   nt_bgc_Fep_in=nt_bgc_Fep,     nt_zbgc_frac_in=nt_zbgc_frac, &
    1798             :            nlt_chl_sw_in=nlt_chl_sw,   nlt_bgc_Sil_in=nlt_bgc_Sil,   nlt_zaero_sw_in=nlt_zaero_sw, &
    1799             :            nlt_bgc_N_in=nlt_bgc_N,     nlt_bgc_Nit_in=nlt_bgc_Nit,   nlt_bgc_Am_in=nlt_bgc_Am,   &
    1800             :            nlt_bgc_DMS_in=nlt_bgc_DMS, nlt_bgc_DMSPp_in=nlt_bgc_DMSPp,                           &
    1801             :            nlt_bgc_DMSPd_in=nlt_bgc_DMSPd,                                                       &
    1802             :            nlt_bgc_DIC_in=nlt_bgc_DIC, nlt_bgc_DOC_in=nlt_bgc_DOC,   nlt_bgc_PON_in=nlt_bgc_PON, &
    1803             :            nlt_bgc_DON_in=nlt_bgc_DON, nlt_bgc_Fed_in=nlt_bgc_Fed,   nlt_bgc_Fep_in=nlt_bgc_Fep, &
    1804             :            nlt_bgc_chl_in=nlt_bgc_chl, nlt_bgc_hum_in=nlt_bgc_hum,   nlt_zaero_in=nlt_zaero,     &
    1805          46 :            bio_index_o_in=bio_index_o, bio_index_in=bio_index)
    1806          46 :       call icepack_warnings_flush(nu_diag)
    1807          46 :       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
    1808           0 :           file=__FILE__, line=__LINE__)
    1809             : 
    1810             :       !-----------------------------------------------------------------
    1811             :       ! final consistency checks
    1812             :       !----------------------------------------------------------------- 
    1813          46 :       if (nbtrcr > icepack_max_nbtrcr) then
    1814           0 :          write (nu_diag,*) ' '
    1815           0 :          write (nu_diag,*) 'nbtrcr > icepack_max_nbtrcr'
    1816           0 :          write (nu_diag,*) 'nbtrcr, icepack_max_nbtrcr:',nbtrcr, icepack_max_nbtrcr
    1817           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1818             :       endif
    1819          46 :       if (.NOT. dEdd_algae) nbtrcr_sw = 1
    1820             : 
    1821          46 :       if (nbtrcr_sw > max_nsw) then
    1822           0 :          write (nu_diag,*) ' '
    1823           0 :          write (nu_diag,*) 'nbtrcr_sw > max_nsw'
    1824           0 :          write (nu_diag,*) 'nbtrcr_sw, max_nsw:',nbtrcr_sw, max_nsw
    1825           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1826             :       endif
    1827             : 
    1828          46 :       if (ntrcr > max_ntrcr) then
    1829           0 :          write(nu_diag,*) 'max_ntrcr < number of namelist tracers'
    1830           0 :          write(nu_diag,*) 'max_ntrcr = ',max_ntrcr,' ntrcr = ',ntrcr
    1831           0 :          call icedrv_system_abort(file=__FILE__,line=__LINE__)
    1832             :       endif                               
    1833             : 
    1834             :       !-----------------------------------------------------------------
    1835             :       ! spew
    1836             :       !-----------------------------------------------------------------
    1837          46 :       if (skl_bgc) then
    1838             : 
    1839           1 :          write(nu_diag,1010) ' skl_bgc                   = ', skl_bgc
    1840           1 :          write(nu_diag,1030) ' bgc_flux_type             = ', bgc_flux_type
    1841           1 :          write(nu_diag,1010) ' restore_bgc               = ', restore_bgc
    1842           1 :          write(nu_diag,*)    ' bgc_data_type             = ', &
    1843           2 :                                trim(bgc_data_type)
    1844           1 :          write(nu_diag,1020) ' number of bio tracers     = ', nbtrcr
    1845           1 :          write(nu_diag,1020) ' number of Isw tracers     = ', nbtrcr_sw
    1846           1 :          write(nu_diag,1020) ' number of autotrophs      = ', n_algae
    1847           1 :          write(nu_diag,1020) ' number of doc          = ', n_doc
    1848           1 :          write(nu_diag,1020) ' number of dic          = ', n_dic
    1849           1 :          write(nu_diag,1020) ' number of don          = ', n_don
    1850           1 :          write(nu_diag,1020) ' number of fed          = ', n_fed
    1851           1 :          write(nu_diag,1020) ' number of fep          = ', n_fep
    1852           1 :          write(nu_diag,1010) ' tr_bgc_N               = ', tr_bgc_N
    1853           1 :          write(nu_diag,1010) ' tr_bgc_C               = ', tr_bgc_C
    1854           1 :          write(nu_diag,1010) ' tr_bgc_chl             = ', tr_bgc_chl
    1855           1 :          write(nu_diag,1010) ' tr_bgc_Nit             = ', tr_bgc_Nit
    1856           1 :          write(nu_diag,1010) ' tr_bgc_Am              = ', tr_bgc_Am
    1857           1 :          write(nu_diag,1010) ' tr_bgc_Sil             = ', tr_bgc_Sil
    1858           1 :          write(nu_diag,1010) ' tr_bgc_hum             = ', tr_bgc_hum
    1859           1 :          write(nu_diag,1010) ' tr_bgc_DMS             = ', tr_bgc_DMS
    1860           1 :          write(nu_diag,1010) ' tr_bgc_PON             = ', tr_bgc_PON
    1861           1 :          write(nu_diag,1010) ' tr_bgc_DON             = ', tr_bgc_DON
    1862           1 :          write(nu_diag,1010) ' tr_bgc_Fe              = ', tr_bgc_Fe 
    1863             :         
    1864          45 :       elseif (z_tracers) then
    1865             : 
    1866           4 :          write(nu_diag,*)    ' bgc_data_type             = ', &
    1867           8 :                                trim(bgc_data_type)
    1868           4 :          write(nu_diag,1010) ' dEdd_algae                = ', dEdd_algae  
    1869           4 :          write(nu_diag,1010) ' modal_aero                = ', modal_aero  
    1870           4 :          write(nu_diag,1010) ' scale_bgc                 = ', scale_bgc
    1871           4 :          write(nu_diag,1010) ' solve_zbgc                = ', solve_zbgc
    1872           4 :          write(nu_diag,1020) ' number of ztracers        = ', nbtrcr
    1873           4 :          write(nu_diag,1020) ' number of Isw tracers     = ', nbtrcr_sw
    1874           4 :          write(nu_diag,1020) ' number of autotrophs      = ', n_algae
    1875           4 :          write(nu_diag,1020) ' number of doc             = ', n_doc
    1876           4 :          write(nu_diag,1020) ' number of dic             = ', n_dic
    1877           4 :          write(nu_diag,1020) ' number of fed             = ', n_fed
    1878           4 :          write(nu_diag,1020) ' number of fep             = ', n_fep
    1879           4 :          write(nu_diag,1020) ' number of aerosols        = ', n_zaero
    1880           4 :          write(nu_diag,1010) ' tr_zaero                  = ', tr_zaero
    1881           4 :          write(nu_diag,1010) ' tr_bgc_Nit                = ', tr_bgc_Nit
    1882           4 :          write(nu_diag,1010) ' tr_bgc_N                  = ', tr_bgc_N
    1883           4 :          write(nu_diag,1010) ' tr_bgc_Am                 = ', tr_bgc_Am
    1884           4 :          write(nu_diag,1010) ' tr_bgc_C                  = ', tr_bgc_C
    1885           4 :          write(nu_diag,1010) ' tr_bgc_Sil                = ', tr_bgc_Sil
    1886           4 :          write(nu_diag,1010) ' tr_bgc_hum                = ', tr_bgc_hum
    1887           4 :          write(nu_diag,1010) ' tr_bgc_chl                = ', tr_bgc_chl
    1888           4 :          write(nu_diag,1010) ' tr_bgc_DMS                = ', tr_bgc_DMS
    1889           4 :          write(nu_diag,1010) ' tr_bgc_PON                = ', tr_bgc_PON
    1890           4 :          write(nu_diag,1010) ' tr_bgc_DON                = ', tr_bgc_DON
    1891           4 :          write(nu_diag,1010) ' tr_bgc_Fe                 = ', tr_bgc_Fe 
    1892           4 :          write(nu_diag,1000) ' grid_o                    = ', grid_o
    1893           4 :          write(nu_diag,1005) ' l_sk                      = ', l_sk
    1894           4 :          write(nu_diag,1000) ' initbio_frac              = ', initbio_frac
    1895           4 :          write(nu_diag,1000) ' frazil_scav               = ', frazil_scav  
    1896             : 
    1897             :       endif  ! skl_bgc or solve_bgc
    1898             : 
    1899             :  1000    format (a30,2x,f9.2)  ! a30 to align formatted, unformatted statements
    1900             :  1005    format (a30,2x,f9.6)  ! float
    1901             :  1010    format (a30,2x,l6)    ! logical
    1902             :  1020    format (a30,2x,i6)    ! integer
    1903             :  1030    format (a30,   a8)    ! character
    1904             : 
    1905          46 :       end subroutine init_zbgc
    1906             : 
    1907             : !=======================================================================
    1908             : 
    1909          92 :       subroutine init_bgc_trcr(nk,              nt_fbri,       &
    1910             :                                nt_bgc,          nlt_bgc,       &
    1911             :                                bgctype,         nt_depend,     &
    1912             :                                ntrcr,           nbtrcr,        &
    1913          92 :                                bgc_tracer_type, trcr_depend,   &
    1914          92 :                                trcr_base,       n_trcr_strata, &
    1915          92 :                                nt_strata,       bio_index)
    1916             : 
    1917             :       integer (kind=int_kind), intent(in) :: &
    1918             :          nk           , & ! counter
    1919             :          nt_depend    , & ! tracer dependency index
    1920             :          nt_fbri
    1921             : 
    1922             :       integer (kind=int_kind), intent(inout) :: &
    1923             :          ntrcr        , & ! number of tracers
    1924             :          nbtrcr       , & ! number of bio tracers
    1925             :          nt_bgc       , & ! tracer index
    1926             :          nlt_bgc          ! bio tracer index
    1927             : 
    1928             :       integer (kind=int_kind), dimension(:), intent(inout) :: &
    1929             :          trcr_depend  , & ! tracer dependencies
    1930             :          n_trcr_strata, & ! number of underlying tracer layers
    1931             :          bio_index        !
    1932             : 
    1933             :       integer (kind=int_kind), dimension(:,:), intent(inout) :: &
    1934             :          nt_strata        ! indices of underlying tracer layers
    1935             : 
    1936             :       real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
    1937             :          trcr_base        ! = 0 or 1 depending on tracer dependency
    1938             :                           ! argument 2:  (1) aice, (2) vice, (3) vsno
    1939             : 
    1940             :       real (kind=dbl_kind), intent(in) :: &
    1941             :          bgctype          ! bio tracer transport type (mobile vs stationary)
    1942             : 
    1943             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
    1944             :          bgc_tracer_type  ! bio tracer transport type array
    1945             : 
    1946             :       ! local variables
    1947             : 
    1948             :       integer (kind=int_kind) :: &
    1949             :          k         , & ! loop index
    1950             :          n_strata  , & ! temporary values
    1951             :          nt_strata1, & ! 
    1952             :          nt_strata2
    1953             : 
    1954             :       real (kind=dbl_kind) :: &
    1955          52 :          trcr_base1, & ! temporary values
    1956          52 :          trcr_base2, &
    1957          52 :          trcr_base3
    1958             : 
    1959             :       character(len=*), parameter :: subname='(init_bgc_trcr)'
    1960             : 
    1961             :       !--------
    1962             : 
    1963          92 :       nt_bgc = ntrcr + 1 
    1964          92 :       nbtrcr = nbtrcr + 1
    1965          92 :       nlt_bgc = nbtrcr
    1966          92 :       bgc_tracer_type(nbtrcr) = bgctype
    1967             :          
    1968          92 :       if (nk > 1) then 
    1969             :          ! include vertical bgc in snow
    1970         228 :          do k = nk, nk+1
    1971         152 :             ntrcr = ntrcr + 1
    1972         152 :             trcr_depend  (nt_bgc + k  ) = 2 ! snow volume
    1973         152 :             trcr_base    (nt_bgc + k,1) = c0
    1974         152 :             trcr_base    (nt_bgc + k,2) = c0
    1975         152 :             trcr_base    (nt_bgc + k,3) = c1
    1976         152 :             n_trcr_strata(nt_bgc + k  ) = 0
    1977         152 :             nt_strata    (nt_bgc + k,1) = 0
    1978         228 :             nt_strata    (nt_bgc + k,2) = 0
    1979             :          enddo
    1980             : 
    1981          76 :          trcr_base1 = c0      
    1982          76 :          trcr_base2 = c1     
    1983          76 :          trcr_base3 = c0
    1984          76 :          n_strata = 1    
    1985          76 :          nt_strata1 = nt_fbri
    1986          76 :          nt_strata2 = 0
    1987             :       else  ! nk = 1
    1988          16 :          trcr_base1 = c1
    1989          16 :          trcr_base2 = c0
    1990          16 :          trcr_base3 = c0
    1991          16 :          n_strata = 0
    1992          16 :          nt_strata1 = 0
    1993          16 :          nt_strata2 = 0
    1994             :       endif ! nk
    1995             : 
    1996         716 :       do k = 1, nk     !in ice
    1997         624 :          ntrcr = ntrcr + 1
    1998         624 :          trcr_depend  (nt_bgc + k - 1  ) = nt_depend
    1999         624 :          trcr_base    (nt_bgc + k - 1,1) = trcr_base1
    2000         624 :          trcr_base    (nt_bgc + k - 1,2) = trcr_base2
    2001         624 :          trcr_base    (nt_bgc + k - 1,3) = trcr_base3
    2002         624 :          n_trcr_strata(nt_bgc + k - 1  ) = n_strata
    2003         624 :          nt_strata    (nt_bgc + k - 1,1) = nt_strata1
    2004         716 :          nt_strata    (nt_bgc + k - 1,2) = nt_strata2
    2005             :       enddo
    2006             : 
    2007          92 :       bio_index (nlt_bgc) = nt_bgc
    2008             : 
    2009          92 :       end subroutine init_bgc_trcr
    2010             : 
    2011             : !=======================================================================
    2012             : 
    2013             :       end module icedrv_init_column
    2014             : 
    2015             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd