LCOV - code coverage report
Current view: top level - columnphysics - icepack_aerosol.F90 (source / functions) Hit Total Coverage
Test: 200804-015008:4c42a82e3d:3:base,travis,quick Lines: 294 436 67.43 %
Date: 2020-08-03 20:10:57 Functions: 2 2 100.00 %

          Line data    Source code
       1             : !=======================================================================
       2             : 
       3             : ! Aerosol tracer within sea ice
       4             : !
       5             : ! authors Marika Holland, NCAR
       6             : !         David Bailey, NCAR
       7             : 
       8             :       module icepack_aerosol
       9             : 
      10             :       use icepack_kinds
      11             :       use icepack_parameters, only: c0, c1, c2, puny, rhoi, rhos, hs_min
      12             :       use icepack_parameters, only: hi_ssl, hs_ssl
      13             :       use icepack_tracers, only: max_aero
      14             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      15             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      16             : 
      17             :       use icepack_zbgc_shared, only: kscavz 
      18             : 
      19             :       implicit none
      20             : 
      21             :       private
      22             :       public :: update_aerosol, update_snow_bgc
      23             : 
      24             : !=======================================================================
      25             : 
      26             :       contains
      27             : 
      28             : !=======================================================================
      29             : 
      30             : !  Increase aerosol in ice or snow surface due to deposition
      31             : !  and vertical cycling
      32             : 
      33      229765 :       subroutine update_aerosol(dt,                   &
      34             :                                 nilyr,    nslyr,      &
      35             :                                 n_aero,     &
      36             :                                 meltt,    melts,      &
      37             :                                 meltb,    congel,     &
      38             :                                 snoice,               &
      39             :                                 fsnow,                &
      40      229765 :                                 aerosno,  aeroice,    &
      41             :                                 aice_old,             &
      42             :                                 vice_old, vsno_old,   &
      43             :                                 vicen, vsnon, aicen,  &
      44      459530 :                                 faero_atm, faero_ocn)
      45             : 
      46             :       integer (kind=int_kind), intent(in) :: &
      47             :          nilyr, nslyr, n_aero
      48             : 
      49             :       real (kind=dbl_kind), intent(in) :: &
      50             :          dt,       & ! time step
      51             :          meltt,    & ! thermodynamic melt/growth rates
      52             :          melts,    &
      53             :          meltb,    &
      54             :          congel,   &
      55             :          snoice,   &
      56             :          fsnow,    &
      57             :          vicen,    & ! ice volume (m)
      58             :          vsnon,    & ! snow volume (m)
      59             :          aicen,    & ! ice area fraction
      60             :          aice_old, & ! values prior to thermodynamic changes
      61             :          vice_old, &
      62             :          vsno_old 
      63             : 
      64             :       real (kind=dbl_kind), dimension(:), &
      65             :          intent(in) :: &
      66             :          faero_atm   ! aerosol deposition rate (W/m^2 s)
      67             : 
      68             :       real (kind=dbl_kind), dimension(:), &
      69             :          intent(inout) :: &
      70             :          faero_ocn   ! aerosol flux to ocean (W/m^2 s)
      71             : 
      72             :       real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
      73             :          aerosno,  aeroice    ! kg/m^2
      74             : 
      75             :       !  local variables
      76             :       integer (kind=int_kind) :: k
      77             : 
      78             :       real (kind=dbl_kind) :: &
      79           0 :          dzssl,  dzssl_new,      & ! snow ssl thickness
      80           0 :          dzint,                  & ! snow interior thickness
      81           0 :          dzssli, dzssli_new,     & ! ice ssl thickness
      82           0 :          dzinti,                 & ! ice interior thickness
      83           0 :          dznew,                  & ! tracks thickness changes
      84           0 :          hs, hi,                 & ! snow/ice thickness (m)
      85           0 :          dhs_evap, dhi_evap,     & ! snow/ice thickness change due to evap
      86           0 :          dhs_melts, dhi_meltt,   & ! ... due to surface melt
      87           0 :          dhs_snoice, dhi_snoice, & ! ... due to snow-ice formation
      88           0 :          dhi_congel, dhi_meltb,  & ! ... due to bottom growth, melt
      89           0 :          hslyr, hilyr,           & ! snow, ice layer thickness (m)
      90           0 :          hslyr_old, hilyr_old,   & ! old snow, ice layer thickness (m)
      91           0 :          hs_old, hi_old,         & ! old snow, ice thickness (m)
      92           0 :          sloss1, sloss2,         & ! aerosol mass loss (kg/m^2)
      93           0 :          ar                        ! 1/aicen(i,j)
      94             : 
      95             :       real (kind=dbl_kind), dimension(max_aero) :: &
      96             :          kscav, kscavsi       ! scavenging by melt water
      97             : 
      98             :       real (kind=dbl_kind), dimension(n_aero) :: &
      99      459530 :          aerotot, aerotot0, & ! for conservation check
     100      459530 :          focn_old             ! for conservation check
     101             : 
     102             :       ! echmod:  this assumes max_aero=6
     103             :       data kscav   / .03_dbl_kind, .20_dbl_kind, .02_dbl_kind, &
     104             :                      .02_dbl_kind, .01_dbl_kind, .01_dbl_kind / 
     105             :       data kscavsi / .03_dbl_kind, .20_dbl_kind, .02_dbl_kind, &
     106             :                      .02_dbl_kind, .01_dbl_kind, .01_dbl_kind / 
     107             : 
     108             :      character(len=*),parameter :: subname='(update_aerosol)'
     109             : 
     110             :     !-------------------------------------------------------------------
     111             :     ! initialize
     112             :     !-------------------------------------------------------------------
     113      459530 :       focn_old(:) = faero_ocn(:)
     114             :       
     115      229765 :       hs_old    = vsno_old/aice_old
     116      229765 :       hi_old    = vice_old/aice_old
     117      229765 :       hslyr_old = hs_old/real(nslyr,kind=dbl_kind)
     118      229765 :       hilyr_old = hi_old/real(nilyr,kind=dbl_kind)
     119             :       
     120      229765 :       dzssl  = min(hslyr_old/c2, hs_ssl)
     121      229765 :       dzssli = min(hilyr_old/c2, hi_ssl)
     122      229765 :       dzint  = hs_old - dzssl
     123      229765 :       dzinti = hi_old - dzssli
     124             :       
     125      229765 :       if (aicen > c0) then
     126      229765 :          ar = c1/aicen
     127             :       else ! ice disappeared during time step
     128           0 :          ar = c1/aice_old
     129             :       endif
     130             : 
     131      229765 :       hs = vsnon*ar
     132      229765 :       hi = vicen*ar
     133      229765 :       dhs_melts  = -melts
     134      229765 :       dhi_snoice = snoice
     135      229765 :       dhs_snoice = dhi_snoice*rhoi/rhos
     136      229765 :       dhi_meltt  = -meltt
     137      229765 :       dhi_meltb  = -meltb
     138      229765 :       dhi_congel = congel
     139             : 
     140             :       dhs_evap = hs - (hs_old + dhs_melts - dhs_snoice &
     141      229765 :                               + fsnow/rhos*dt)
     142             :       dhi_evap = hi - (hi_old + dhi_meltt + dhi_meltb &
     143      229765 :                               + dhi_congel + dhi_snoice)
     144             : 
     145      459530 :       do k = 1, n_aero
     146           0 :          aerotot0(k) = aerosno(k,2) + aerosno(k,1) &
     147      459530 :                      + aeroice(k,2) + aeroice(k,1)
     148             :       enddo
     149             :       
     150             :     !-------------------------------------------------------------------
     151             :     ! evaporation
     152             :     !-------------------------------------------------------------------
     153      229765 :       dzint  = dzint  + min(dzssl  + dhs_evap, c0)
     154      229765 :       dzinti = dzinti + min(dzssli + dhi_evap, c0)
     155      229765 :       dzssl  = max(dzssl  + dhs_evap, c0)
     156      229765 :       dzssli = max(dzssli + dhi_evap, c0)
     157             : 
     158             :     !-------------------------------------------------------------------
     159             :     ! basal ice growth
     160             :     !-------------------------------------------------------------------
     161      229765 :       dzinti = dzinti + dhi_congel
     162             : 
     163             :     !-------------------------------------------------------------------
     164             :     ! surface snow melt
     165             :     !-------------------------------------------------------------------
     166      229765 :       if (-dhs_melts > puny) then
     167        8104 :          do k = 1, n_aero
     168        4052 :             sloss1 = c0
     169        4052 :             sloss2 = c0
     170        4052 :             if (dzssl > puny)  &
     171           0 :                  sloss1 = kscav(k)*aerosno(k,1)  &
     172        3687 :                                   *min(-dhs_melts,dzssl)/dzssl
     173        4052 :             aerosno(k,1) = aerosno(k,1) - sloss1
     174        4052 :             if (dzint > puny)  &
     175           0 :                  sloss2 = kscav(k)*aerosno(k,2) &
     176        4052 :                                   *max(-dhs_melts-dzssl,c0)/dzint
     177        4052 :             aerosno(k,2) = aerosno(k,2) - sloss2
     178        8104 :             faero_ocn(k) = faero_ocn(k) + (sloss1+sloss2)/dt
     179             :          enddo  ! n_aero
     180             : 
     181             :          ! update snow thickness
     182        4052 :          dzint=dzint+min(dzssl+dhs_melts, c0)
     183        4052 :          dzssl=max(dzssl+dhs_melts, c0)
     184             : 
     185        4052 :          if ( dzssl <= puny ) then ! ssl melts away
     186         768 :             aerosno(:,2) = aerosno(:,1) + aerosno(:,2)
     187         768 :             aerosno(:,1) = c0
     188         384 :             dzssl = max(dzssl, c0)
     189             :          endif
     190        4052 :          if (dzint <= puny ) then  ! all snow melts away
     191           0 :             aeroice(:,1) = aeroice(:,1) &
     192         752 :                          + aerosno(:,1) + aerosno(:,2)
     193        1880 :             aerosno(:,:) = c0
     194         376 :             dzint = max(dzint, c0)
     195             :          endif
     196             :       endif
     197             : 
     198             :     !-------------------------------------------------------------------
     199             :     ! surface ice melt
     200             :     !-------------------------------------------------------------------
     201      229765 :       if (-dhi_meltt > puny) then
     202       18234 :          do k = 1, n_aero
     203        9117 :             sloss1 = c0
     204        9117 :             sloss2 = c0
     205        9117 :             if (dzssli > puny)  &
     206           0 :                  sloss1 = kscav(k)*aeroice(k,1)  &
     207        9117 :                                   *min(-dhi_meltt,dzssli)/dzssli
     208        9117 :             aeroice(k,1) = aeroice(k,1) - sloss1
     209        9117 :             if (dzinti > puny)  &
     210           0 :                  sloss2 = kscav(k)*aeroice(k,2)  &
     211        9117 :                                   *max(-dhi_meltt-dzssli,c0)/dzinti
     212        9117 :             aeroice(k,2) = aeroice(k,2) - sloss2
     213       18234 :             faero_ocn(k) = faero_ocn(k) + (sloss1+sloss2)/dt
     214             :          enddo
     215             :          
     216        9117 :          dzinti = dzinti + min(dzssli+dhi_meltt, c0)
     217        9117 :          dzssli = max(dzssli+dhi_meltt, c0)
     218        9117 :          if (dzssli <= puny) then   ! ssl ice melts away
     219          30 :             do k = 1, n_aero
     220          15 :                aeroice(k,2) = aeroice(k,1) + aeroice(k,2)
     221          30 :                aeroice(k,1) = c0
     222             :             enddo
     223          15 :             dzssli = max(dzssli, c0)
     224             :          endif
     225        9117 :          if (dzinti <= puny) then   ! all ice melts away
     226           0 :             do k = 1, n_aero
     227           0 :                faero_ocn(k) = faero_ocn(k)  &
     228           0 :                             + (aeroice(k,1)+aeroice(k,2))/dt
     229           0 :                aeroice(k,:)=c0
     230             :             enddo
     231           0 :             dzinti = max(dzinti, c0)
     232             :          endif
     233             :       endif
     234             : 
     235             :     !-------------------------------------------------------------------
     236             :     ! basal ice melt.  Assume all aero lost in basal melt
     237             :     !-------------------------------------------------------------------
     238      229765 :       if (-dhi_meltb > puny) then
     239      217908 :          do k=1,n_aero
     240      108954 :             sloss1=c0
     241      108954 :             sloss2=c0
     242      108954 :             if (dzssli > puny)  &
     243             :                  sloss1 = max(-dhi_meltb-dzinti, c0)  &
     244      108939 :                         *aeroice(k,1)/dzssli
     245      108954 :             aeroice(k,1) = aeroice(k,1) - sloss1
     246      108954 :             if (dzinti > puny)  &
     247             :                  sloss2 = min(-dhi_meltb, dzinti)  &
     248      108954 :                         *aeroice(k,2)/dzinti
     249      108954 :             aeroice(k,2) = aeroice(k,2) - sloss2
     250      217908 :             faero_ocn(k) = faero_ocn(k) + (sloss1+sloss2)/dt
     251             :          enddo
     252             : 
     253      108954 :          dzssli = dzssli + min(dzinti+dhi_meltb, c0)
     254      108954 :          dzinti = max(dzinti+dhi_meltb, c0)           
     255             :       endif
     256             : 
     257             :     !-------------------------------------------------------------------
     258             :     ! snowfall
     259             :     !-------------------------------------------------------------------
     260      229765 :       if (fsnow > c0) dzssl = dzssl + fsnow/rhos*dt
     261             : 
     262             :     !-------------------------------------------------------------------
     263             :     ! snow-ice formation
     264             :     !-------------------------------------------------------------------
     265      229765 :       if (dhs_snoice > puny) then
     266        6600 :          do k = 1, n_aero
     267        3300 :             sloss1 = c0
     268        3300 :             sloss2 = c0
     269        3300 :             if (dzint > puny)  &
     270             :                  sloss2 = min(dhs_snoice, dzint)  &
     271        3300 :                         *aerosno(k,2)/dzint
     272        3300 :             aerosno(k,2) = aerosno(k,2) - sloss2
     273        3300 :             if (dzssl > puny)  &
     274             :                  sloss1 = max(dhs_snoice-dzint, c0)  &
     275        3300 :                         *aerosno(k,1)/dzssl
     276        3300 :             aerosno(k,1) = aerosno(k,1) - sloss1
     277           0 :             aeroice(k,1) = aeroice(k,1) &
     278        3300 :                          + (c1-kscavsi(k))*(sloss2+sloss1)
     279           0 :             faero_ocn(k) = faero_ocn(k) &
     280        6600 :                          + kscavsi(k)*(sloss2+sloss1)/dt
     281             :          enddo
     282        3300 :          dzssl  = dzssl - max(dhs_snoice-dzint, c0)
     283        3300 :          dzint  = max(dzint-dhs_snoice, c0)
     284        3300 :          dzssli = dzssli + dhi_snoice
     285             :       endif
     286             : 
     287             :     !-------------------------------------------------------------------
     288             :     ! aerosol deposition
     289             :     !-------------------------------------------------------------------
     290      229765 :       if (aicen > c0) then
     291      229765 :          hs = vsnon * ar
     292             :       else
     293           0 :          hs = c0
     294             :       endif
     295      229765 :       if (hs > hs_min) then    ! should this really be hs_min or 0? 
     296             :          ! should use same hs_min value as in radiation
     297      378830 :          do k=1,n_aero
     298           0 :             aerosno(k,1) = aerosno(k,1) &
     299      378830 :                          + faero_atm(k)*dt*aicen
     300             :          enddo
     301             :       else
     302       80700 :          do k=1,n_aero
     303           0 :             aeroice(k,1) = aeroice(k,1) &
     304       80700 :                          + faero_atm(k)*dt*aicen
     305             :          enddo
     306             :       endif
     307             : 
     308             :     !-------------------------------------------------------------------
     309             :     ! redistribute aerosol within vertical layers
     310             :     !-------------------------------------------------------------------
     311      229765 :       if (aicen > c0) then
     312      229765 :          hs = vsnon  * ar     ! new snow thickness
     313      229765 :          hi = vicen  * ar     ! new ice thickness
     314             :       else
     315           0 :          hs = c0
     316           0 :          hi = c0
     317             :       endif
     318      229765 :       if (dzssl <= puny) then   ! nothing in SSL
     319        6372 :          do k=1,n_aero
     320        3186 :             aerosno(k,2) = aerosno(k,2) + aerosno(k,1)
     321        6372 :             aerosno(k,1) = c0
     322             :          enddo
     323             :       endif
     324      229765 :       if (dzint <= puny) then   ! nothing in Snow Int
     325       78172 :          do k = 1, n_aero
     326       39086 :             aeroice(k,1) = aeroice(k,1) + aerosno(k,2)
     327       78172 :             aerosno(k,2) = c0
     328             :          enddo
     329             :       endif
     330      229765 :       if (dzssli <= puny) then  ! nothing in Ice SSL
     331          30 :          do k = 1, n_aero
     332          15 :             aeroice(k,2) = aeroice(k,2) + aeroice(k,1)
     333          30 :             aeroice(k,1) = c0
     334             :          enddo
     335             :       endif
     336             :       
     337      229765 :       if (dzinti <= puny) then  ! nothing in Ice INT
     338           0 :          do k = 1, n_aero
     339           0 :             faero_ocn(k) = faero_ocn(k) &
     340           0 :                          + (aeroice(k,1)+aeroice(k,2))/dt
     341           0 :             aeroice(k,:)=c0
     342             :          enddo
     343             :       endif
     344             :       
     345      229765 :       hslyr      = hs/real(nslyr,kind=dbl_kind)
     346      229765 :       hilyr      = hi/real(nilyr,kind=dbl_kind)
     347      229765 :       dzssl_new  = min(hslyr/c2, hs_ssl)
     348      229765 :       dzssli_new = min(hilyr/c2, hi_ssl)
     349             : 
     350      229765 :       if (hs > hs_min) then
     351      378830 :          do k = 1, n_aero
     352      189415 :             dznew = min(dzssl_new-dzssl, c0)
     353      189415 :             sloss1 = c0
     354      189415 :             if (dzssl > puny) &
     355      189415 :                  sloss1 = dznew*aerosno(k,1)/dzssl ! not neccesarily a loss
     356      189415 :             dznew = max(dzssl_new-dzssl, c0)
     357      189415 :             if (dzint > puny) &
     358      189411 :                  sloss1 = sloss1 + aerosno(k,2)*dznew/dzint
     359      189415 :             aerosno(k,1) = aerosno(k,1) + sloss1 
     360      378830 :             aerosno(k,2) = aerosno(k,2) - sloss1
     361             :          enddo
     362             :       else
     363           0 :          aeroice(:,1) = aeroice(:,1)  &
     364       80700 :                       + aerosno(:,1) + aerosno(:,2)
     365      201750 :          aerosno(:,:) = c0
     366             :       endif
     367             :       
     368      229765 :       if (vicen > puny) then ! may want a limit on hi instead?
     369      459478 :          do k = 1, n_aero
     370      229739 :             sloss2 = c0
     371      229739 :             dznew = min(dzssli_new-dzssli, c0)
     372      229739 :             if (dzssli > puny) & 
     373      229724 :                  sloss2 = dznew*aeroice(k,1)/dzssli
     374      229739 :             dznew = max(dzssli_new-dzssli, c0)
     375      229739 :             if (dzinti > puny) & 
     376      229739 :                  sloss2 = sloss2 + aeroice(k,2)*dznew/dzinti
     377      229739 :             aeroice(k,1) = aeroice(k,1) + sloss2 
     378      459478 :             aeroice(k,2) = aeroice(k,2) - sloss2
     379             :          enddo
     380             :       else
     381          52 :          faero_ocn(:) = faero_ocn(:) + (aeroice(:,1)+aeroice(:,2))/dt
     382         130 :          aeroice(:,:) = c0
     383             :       endif
     384             :       
     385             :     !-------------------------------------------------------------------
     386             :     ! check conservation
     387             :     !-------------------------------------------------------------------
     388      459530 :       do k = 1, n_aero
     389           0 :          aerotot(k) = aerosno(k,2) + aerosno(k,1) &
     390      229765 :                     + aeroice(k,2) + aeroice(k,1)
     391      229765 :          if ((aerotot(k)-aerotot0(k)) &
     392           0 :               - (   faero_atm(k)*aicen &
     393      229765 :               - (faero_ocn(k)-focn_old(k)) )*dt  > puny) then
     394             : 
     395           0 :             write(warnstr,*) subname, 'aerosol tracer:  ',k
     396           0 :             call icepack_warnings_add(warnstr)
     397           0 :             write(warnstr,*) subname, 'aerotot-aerotot0 ',aerotot(k)-aerotot0(k)
     398           0 :             call icepack_warnings_add(warnstr)
     399           0 :             write(warnstr,*) subname, 'faero_atm-faero_ocn      ', &
     400           0 :                  (faero_atm(k)*aicen-(faero_ocn(k)-focn_old(k)))*dt
     401           0 :             call icepack_warnings_add(warnstr)
     402             :          endif
     403             :       enddo
     404             : 
     405             :     !-------------------------------------------------------------------
     406             :     ! check for negative values
     407             :     !-------------------------------------------------------------------
     408             : 
     409             : !echmod:  note that this does not test or fix all aero tracers         
     410           0 :       if (aeroice(1,1) < -puny .or. &
     411           0 :           aeroice(1,2) < -puny .or. &
     412      229765 :           aerosno(1,1) < -puny .or. &
     413             :           aerosno(1,2) < -puny) then
     414             : 
     415           0 :          write(warnstr,*) subname, 'aerosol negative in aerosol code'
     416           0 :          call icepack_warnings_add(warnstr)
     417             : 
     418           0 :          aeroice(1,1) = max(aeroice(1,1), c0)
     419           0 :          aeroice(1,2) = max(aeroice(1,2), c0)
     420           0 :          aerosno(1,1) = max(aerosno(1,1), c0)
     421           0 :          aerosno(1,2) = max(aerosno(1,2), c0)
     422             : 
     423             :       endif
     424             : 
     425      229765 :       end subroutine update_aerosol
     426             : 
     427             : !=======================================================================
     428             : 
     429             : !  Increase aerosol in snow surface due to deposition
     430             : !  and vertical cycling : after update_aerosol
     431             : 
     432      199620 :       subroutine update_snow_bgc (dt,     nblyr,       &
     433             :                                 nslyr,                 &
     434             :                                 meltt,    melts,       &
     435             :                                 meltb,    congel,      &
     436             :                                 snoice,   nbtrcr,      &
     437             :                                 fsnow,    ntrcr,       &
     438      199620 :                                 trcrn,    bio_index,   &
     439      199620 :                                 aice_old, zbgc_snow,   &
     440             :                                 vice_old, vsno_old,    &
     441             :                                 vicen,    vsnon,       &
     442      199620 :                                 aicen,    flux_bio_atm,&
     443      199620 :                                 zbgc_atm, flux_bio)
     444             : 
     445             :       integer (kind=int_kind), intent(in) :: &
     446             :          nbtrcr,             & ! number of distinct snow tracers
     447             :          nblyr,              & ! number of bio layers
     448             :          nslyr,              & ! number of snow layers
     449             :          ntrcr                 ! number of tracers
     450             : 
     451             :       integer (kind=int_kind), dimension (nbtrcr), intent(in) :: &
     452             :          bio_index       
     453             : 
     454             :       real (kind=dbl_kind), intent(in) :: &
     455             :          dt                    ! time step
     456             : 
     457             :       real (kind=dbl_kind), intent(in) :: &
     458             :          meltt,    & ! thermodynamic melt/growth rates
     459             :          melts,    &
     460             :          meltb,    &
     461             :          congel,   &
     462             :          snoice,   &
     463             :          fsnow,    &
     464             :          vicen,    & ! ice volume (m)
     465             :          vsnon,    & ! snow volume (m)
     466             :          aicen,    & ! ice area fraction
     467             :          aice_old, & ! values prior to thermodynamic changes
     468             :          vice_old, &
     469             :          vsno_old 
     470             : 
     471             :       real (kind=dbl_kind),dimension(nbtrcr), intent(inout) :: &
     472             :          zbgc_snow, & ! aerosol contribution from snow to ice
     473             :          zbgc_atm,  & ! and atm to ice concentration * volume (kg or mmol/m^3*m)
     474             :          flux_bio     ! total ocean tracer flux (mmol/m^2/s)
     475             : 
     476             :       real (kind=dbl_kind), dimension(nbtrcr), &
     477             :          intent(in) :: &
     478             :          flux_bio_atm   ! aerosol deposition rate (kg or mmol/m^2 s)
     479             : 
     480             :       real (kind=dbl_kind), dimension(ntrcr), &
     481             :          intent(inout) :: &
     482             :          trcrn       ! ice/snow tracer array
     483             : 
     484             :       !  local variables
     485             : 
     486             :       integer (kind=int_kind) ::  k, n
     487             : 
     488             :       real (kind=dbl_kind) :: &
     489       48078 :          dzssl,  dzssl_new,      & ! snow ssl thickness
     490       48078 :          dzint,  dzint_new,      & ! snow interior thickness
     491       48078 :          hs,                     & ! snow thickness (m)
     492       48078 :          dhs_evap,               & ! snow thickness change due to evap
     493       48078 :          dhs_melts,              & ! ... due to surface melt
     494       48078 :          dhs_snoice,             & ! ... due to snow-ice formation
     495       48078 :          hslyr,                  & ! snow layer thickness (m)
     496       48078 :          hslyr_old,              & ! old snow layer thickness (m)
     497       48078 :          hs_old,                 & ! old snow thickness (m)
     498       48078 :          dznew,                  & ! change in the snow sl (m)
     499       48078 :          sloss1, sloss2,         & ! aerosol mass loss (kg/m^2)
     500       48078 :          ar                        ! 1/aicen(i,j)
     501             : 
     502             :       real (kind=dbl_kind), dimension(nbtrcr) :: &
     503     2208194 :          aerotot, aerotot0, & ! for conservation check (mmol/m^3)
     504     1279678 :          aero_cons        , & ! for conservation check (mmol/m^2)
     505     1279678 :          flux_bio_o           ! initial ocean tracer flux (mmol/m^2/s)
     506             : 
     507             :       real (kind=dbl_kind), dimension(nbtrcr,2) :: &
     508     2304350 :          aerosno,  & ! kg/m^2
     509     2352428 :          aerosno0    ! for diagnostic prints
     510             : 
     511             :       character(len=*),parameter :: subname='(update_snow_bgc)'
     512             : 
     513             :     !-------------------------------------------------------------------
     514             :     ! initialize
     515             :     !-------------------------------------------------------------------
     516     8517572 :          aerosno (:,:) = c0
     517     8517572 :          aerosno0(:,:) = c0
     518     4158976 :          aero_cons(:) = c0
     519     4158976 :          zbgc_snow(:) = c0
     520     4158976 :          zbgc_atm(:) = c0
     521             : 
     522      199620 :          hs_old    = vsno_old/aice_old
     523      199620 :          hslyr_old = hs_old/real(nslyr,kind=dbl_kind)
     524             : 
     525      199620 :          dzssl  = min(hslyr_old/c2, hs_ssl)
     526      199620 :          dzint  = hs_old - dzssl
     527             : 
     528      199620 :          if (aicen > c0) then
     529      199620 :             ar = c1/aicen
     530      199620 :             hs = vsnon*ar
     531      199620 :             dhs_melts  = -melts
     532      199620 :             dhs_snoice = snoice*rhoi/rhos
     533             :          else ! ice disappeared during time step
     534           0 :             ar = c1
     535           0 :             hs = vsnon/aice_old
     536           0 :             dhs_melts  = -melts
     537           0 :             dhs_snoice = snoice*rhoi/rhos
     538             :          endif
     539             : 
     540             :          dhs_evap = hs - (hs_old + dhs_melts - dhs_snoice &
     541      199620 :                                  + fsnow/rhos*dt)
     542             : 
     543             :          ! trcrn() has units kg/m^3
     544             : 
     545      199620 :       if ((vsno_old .le. puny) .or. (vsnon .le. puny)) then
     546             : 
     547       35012 :          do k=1,nbtrcr
     548       30736 :             flux_bio(k) = flux_bio(k) +  &
     549       61472 :                          (trcrn(bio_index(k)+ nblyr+1)*dzssl+ &
     550       61872 :                           trcrn(bio_index(k)+ nblyr+2)*dzint)/dt  
     551       31136 :             trcrn(bio_index(k) + nblyr+1) = c0
     552       31136 :             trcrn(bio_index(k) + nblyr+2) = c0
     553       30736 :             zbgc_atm(k) = zbgc_atm(k) &
     554       33074 :                                 + flux_bio_atm(k)*dt 
     555             :          enddo
     556             : 
     557             :       else 
     558             :          
     559     4125902 :          do k=1,nbtrcr
     560     3928220 :             flux_bio_o(k) = flux_bio(k)
     561     3928220 :             aerosno (k,1) = trcrn(bio_index(k)+ nblyr+1) * dzssl
     562     3928220 :             aerosno (k,2) = trcrn(bio_index(k)+ nblyr+2) * dzint
     563    11784660 :             aerosno0(k,:) = aerosno(k,:)
     564     4125902 :             aerotot0(k)   = aerosno(k,2) + aerosno(k,1) 
     565             :          enddo
     566             : 
     567             :     !-------------------------------------------------------------------
     568             :     ! evaporation
     569             :     !-------------------------------------------------------------------
     570      197682 :          dzint  = dzint  + min(dzssl  + dhs_evap, c0)
     571      197682 :          dzssl  = max(dzssl  + dhs_evap, c0)
     572             : 
     573             :     !-------------------------------------------------------------------
     574             :     ! surface snow melt
     575             :     !-------------------------------------------------------------------
     576      197682 :          if (-dhs_melts > puny) then
     577      202290 :             do k = 1, nbtrcr
     578      192576 :                sloss1 = c0
     579      192576 :                sloss2 = c0
     580      192576 :                if (dzssl > puny)  &
     581       54096 :                   sloss1 = kscavz(k)*aerosno(k,1)  &
     582      192576 :                                    *min(-dhs_melts,dzssl)/dzssl
     583      192576 :                aerosno(k,1) = aerosno(k,1) - sloss1
     584      192576 :                if (dzint > puny)  &
     585       54096 :                    sloss2 = kscavz(k)*aerosno(k,2) &
     586      192576 :                                     *max(-dhs_melts-dzssl,c0)/dzint
     587      192576 :                aerosno(k,2) = aerosno(k,2) - sloss2
     588      202290 :                zbgc_snow(k) = zbgc_snow(k) + (sloss1+sloss2)
     589             :             enddo  ! 
     590             : 
     591             :             ! update snow thickness
     592        9714 :             dzint=dzint+min(dzssl+dhs_melts, c0)
     593        9714 :             dzssl=max(dzssl+dhs_melts, c0)
     594             : 
     595        9714 :             if ( dzssl <= puny ) then ! ssl melts away
     596           0 :                aerosno(:,2) = aerosno(:,1) + aerosno(:,2)
     597           0 :                aerosno(:,1) = c0
     598           0 :                dzssl = max(dzssl, c0)
     599             :             endif
     600        9714 :             if (dzint <= puny ) then  ! all snow melts away
     601           0 :                zbgc_snow(:) = zbgc_snow(:) &
     602           0 :                                 + max(c0,aerosno(:,1) + aerosno(:,2))
     603           0 :                aerosno(:,:) = c0
     604           0 :                dzint = max(dzint, c0)
     605             :             endif
     606             :          endif
     607             : 
     608             :     !-------------------------------------------------------------------
     609             :     ! snowfall
     610             :     !-------------------------------------------------------------------
     611      197682 :          if (fsnow > c0) dzssl = dzssl + fsnow/rhos*dt 
     612             : 
     613             :     !-------------------------------------------------------------------
     614             :     ! snow-ice formation
     615             :     !-------------------------------------------------------------------
     616      197682 :          if (dhs_snoice > puny) then
     617      392826 :             do k = 1, nbtrcr
     618      374120 :                sloss1 = c0
     619      374120 :                sloss2 = c0
     620      374120 :                if (dzint > puny)  &
     621             :                   sloss2 = min(dhs_snoice, dzint)  &
     622      374120 :                            *aerosno(k,2)/dzint
     623      374120 :                aerosno(k,2) = aerosno(k,2) - sloss2
     624      374120 :                if (dzssl > puny)  &
     625             :                   sloss1 = max(dhs_snoice-dzint, c0)  &
     626      374120 :                            *aerosno(k,1)/dzssl
     627      374120 :                aerosno(k,1) = aerosno(k,1) - sloss1
     628      112740 :                zbgc_snow(k) = zbgc_snow(k) &
     629      392826 :                                + (sloss2+sloss1)
     630             :             enddo
     631       18706 :             dzssl  = dzssl - max(dhs_snoice-dzint, c0)
     632       18706 :             dzint  = max(dzint-dhs_snoice, c0)
     633             :          endif
     634             : 
     635             :     !-------------------------------------------------------------------
     636             :     ! aerosol deposition
     637             :     !-------------------------------------------------------------------
     638      197682 :          if (aicen > c0) then
     639      197682 :             hs = vsnon * ar
     640             :          else
     641           0 :             hs = c0
     642             :          endif
     643      197682 :          if (hs >= hs_min)  then !should this really be hs_min or 0? 
     644             :                                   ! should use same hs_min value as in radiation
     645     4090242 :             do k=1,nbtrcr
     646      877680 :                aerosno(k,1) = aerosno(k,1) &
     647     4090242 :                                 + flux_bio_atm(k)*dt
     648             :             enddo
     649             :          else  
     650       35660 :             do k=1,nbtrcr
     651       20100 :                zbgc_atm(k) = zbgc_atm(k) &
     652       35660 :                                 + flux_bio_atm(k)*dt
     653             :             enddo
     654             :          endif
     655             : 
     656             :     !-------------------------------------------------------------------
     657             :     ! redistribute aerosol within vertical layers
     658             :     !-------------------------------------------------------------------
     659      197682 :          if (aicen > c0) then
     660      197682 :             hs = vsnon  * ar     ! new snow thickness
     661             :          else
     662           0 :             hs = c0
     663             :          endif
     664      197682 :          if (dzssl <= puny) then   ! nothing in SSL
     665           0 :             do k=1,nbtrcr
     666           0 :                aerosno(k,2) = aerosno(k,2) + aerosno(k,1)
     667           0 :                aerosno(k,1) = c0
     668             :             enddo
     669             :          endif
     670      197682 :          if (dzint <= puny) then   ! nothing in Snow Int
     671       67000 :             do k = 1, nbtrcr
     672       63652 :                zbgc_snow(k) = zbgc_snow(k) + max(c0,aerosno(k,2))
     673       67000 :                aerosno(k,2) = c0
     674             :             enddo
     675             :          endif
     676             : 
     677      197682 :          hslyr      = hs/real(nslyr,kind=dbl_kind)
     678      197682 :          dzssl_new  = min(hslyr/c2, hs_ssl)
     679      197682 :          dzint_new  = hs - dzssl_new
     680             : 
     681      197682 :          if (hs > hs_min) then !should this really be hs_min or 0? 
     682     4090242 :             do k = 1, nbtrcr
     683     3894420 :                dznew = min(dzssl_new-dzssl, c0)
     684     3894420 :                sloss1 = c0
     685     3894420 :                if (dzssl > puny) &
     686     3894420 :                   sloss1 = dznew*aerosno(k,1)/dzssl ! not neccesarily a loss
     687     3894420 :                   dznew = max(dzssl_new-dzssl, c0)
     688     3894420 :                if (dzint > puny) &
     689     3862120 :                   sloss1 = sloss1 + aerosno(k,2)*dznew/dzint
     690     3894420 :                aerosno(k,1) = aerosno(k,1) + sloss1 
     691     4090242 :                aerosno(k,2) = aerosno(k,2) - sloss1
     692             :             enddo
     693             :          else
     694           0 :             zbgc_snow(:) = zbgc_snow(:)  &
     695       35660 :                              + max(c0,aerosno(:,1) + aerosno(:,2))
     696       73180 :             aerosno(:,:) = c0
     697             :          endif
     698             : 
     699             :     !-------------------------------------------------------------------
     700             :     ! check conservation
     701             :     !-------------------------------------------------------------------
     702     4125902 :          do k = 1, nbtrcr
     703      897780 :             aerotot(k) = aerosno(k,2) + aerosno(k,1) &
     704     3928220 :                        + zbgc_snow(k) + zbgc_atm(k) 
     705      897780 :             aero_cons(k) = aerotot(k)-aerotot0(k) &
     706      897780 :                           - (    flux_bio_atm(k) & 
     707     3928220 :                           - (flux_bio(k)-flux_bio_o(k))) * dt
     708     4125902 :             if (aero_cons(k)  > puny .or. zbgc_snow(k) + zbgc_atm(k) < c0) then             
     709           0 :                write(warnstr,*) subname, 'Conservation failure: aerosols in snow'
     710           0 :                call icepack_warnings_add(warnstr)
     711           0 :                write(warnstr,*) subname, 'test aerosol 1'
     712           0 :                call icepack_warnings_add(warnstr)
     713           0 :                write(warnstr,*) subname, 'aerosol tracer:  ',k
     714           0 :                call icepack_warnings_add(warnstr)
     715           0 :                write(warnstr,*) subname, 'aero_cons(k),puny:', aero_cons(k),puny
     716           0 :                call icepack_warnings_add(warnstr)
     717           0 :                write(warnstr,*) subname, 'aerotot,aerotot0 ',aerotot(k),aerotot0(k)
     718           0 :                call icepack_warnings_add(warnstr)
     719           0 :                write(warnstr,*) subname, ' aerosno(k,2),aerosno(k,1) ', aerosno(k,2),aerosno(k,1)
     720           0 :                call icepack_warnings_add(warnstr)
     721           0 :                write(warnstr,*) subname, 'flux_bio_atm(k)*aicen*dt', &
     722           0 :                     flux_bio_atm(k)*aicen*dt
     723           0 :                call icepack_warnings_add(warnstr)
     724           0 :                write(warnstr,*) subname, 'zbgc_snow(k)', &
     725           0 :                     zbgc_snow(k)
     726           0 :                call icepack_warnings_add(warnstr)
     727           0 :                write(warnstr,*) subname, 'zbgc_atm(k)', &
     728           0 :                     zbgc_atm(k)
     729           0 :                call icepack_warnings_add(warnstr)
     730             :             endif
     731             :          enddo
     732             : 
     733             :     !-------------------------------------------------------------------
     734             :     ! reload tracers
     735             :     !-------------------------------------------------------------------
     736      197682 :          if (vsnon > puny) then
     737     4125902 :          do k = 1,nbtrcr
     738     3928220 :              trcrn(bio_index(k)+nblyr+1)=aerosno(k,1)/dzssl_new 
     739     4125902 :              trcrn(bio_index(k)+nblyr+2)=aerosno(k,2)/dzint_new
     740             :          enddo
     741             :          else
     742           0 :          do k = 1,nbtrcr
     743           0 :             zbgc_snow(k) = (zbgc_snow(k) + aerosno(k,1) + aerosno(k,2))
     744           0 :             trcrn(bio_index(k)+nblyr+1)= c0
     745           0 :             trcrn(bio_index(k)+nblyr+2)= c0
     746             :          enddo
     747             :          endif
     748             :     !-------------------------------------------------------------------
     749             :     ! check for negative values
     750             :     !-------------------------------------------------------------------
     751     8100282 :          if (minval(aerosno(:,1)) < -puny  .or. &
     752       46160 :             minval(aerosno(:,2)) < -puny) then
     753             : 
     754           0 :             write(warnstr,*) subname, 'Snow aerosol negative in update_snow_bgc'
     755           0 :             call icepack_warnings_add(warnstr)
     756           0 :             write(warnstr,*) subname, 'aicen= '      ,aicen
     757           0 :             call icepack_warnings_add(warnstr)
     758           0 :             write(warnstr,*) subname, 'vicen= '      ,vicen
     759           0 :             call icepack_warnings_add(warnstr)
     760           0 :             write(warnstr,*) subname, 'vsnon= '      ,vsnon
     761           0 :             call icepack_warnings_add(warnstr)
     762           0 :             write(warnstr,*) subname, 'viceold= '    ,vice_old
     763           0 :             call icepack_warnings_add(warnstr)
     764           0 :             write(warnstr,*) subname, 'vsnoold= '    ,vsno_old
     765           0 :             call icepack_warnings_add(warnstr)
     766           0 :             write(warnstr,*) subname, 'melts= '      ,melts
     767           0 :             call icepack_warnings_add(warnstr)
     768           0 :             write(warnstr,*) subname, 'meltt= '      ,meltt
     769           0 :             call icepack_warnings_add(warnstr)
     770           0 :             write(warnstr,*) subname, 'meltb= '      ,meltb
     771           0 :             call icepack_warnings_add(warnstr)
     772           0 :             write(warnstr,*) subname, 'congel= '     ,congel
     773           0 :             call icepack_warnings_add(warnstr)
     774           0 :             write(warnstr,*) subname, 'snoice= '     ,snoice
     775           0 :             call icepack_warnings_add(warnstr)
     776           0 :             write(warnstr,*) subname, 'aero evap snow= '  ,dhs_evap
     777           0 :             call icepack_warnings_add(warnstr)
     778           0 :             write(warnstr,*) subname, 'fsnow= '      ,fsnow
     779           0 :             call icepack_warnings_add(warnstr)
     780           0 :             do k = 1, nbtrcr
     781           0 :               write(warnstr,*) subname, 'NBTRCR value k = ', k
     782           0 :               call icepack_warnings_add(warnstr)
     783           0 :               write(warnstr,*) subname, 'aero snowssl (k)= '    ,aerosno0(k,1)
     784           0 :               call icepack_warnings_add(warnstr)
     785           0 :               write(warnstr,*) subname, 'aero new snowssl (k)= ',aerosno(k,1)
     786           0 :               call icepack_warnings_add(warnstr)
     787           0 :               write(warnstr,*) subname, 'aero snowint (k)= '    ,aerosno0(k,2)
     788           0 :               call icepack_warnings_add(warnstr)
     789           0 :               write(warnstr,*) subname, 'aero new snowint(k)= ',aerosno(k,2)
     790           0 :               call icepack_warnings_add(warnstr)
     791           0 :               write(warnstr,*) subname, 'flux_bio_atm(k)= ' , flux_bio_atm(k)
     792           0 :               call icepack_warnings_add(warnstr)
     793           0 :               write(warnstr,*) subname, 'zbgc_snow(k)= '  ,zbgc_snow(k)
     794           0 :               call icepack_warnings_add(warnstr)
     795           0 :               write(warnstr,*) subname, 'zbgc_atm(k)= '  ,zbgc_atm(k)
     796           0 :               call icepack_warnings_add(warnstr)
     797             : 
     798      197682 :               do n = 1,2
     799           0 :                 trcrn(bio_index(k)+nblyr+n)=max(trcrn(bio_index(k)+nblyr+n), c0)
     800             :               enddo
     801             :             enddo
     802             :          endif
     803             :         endif
     804             : 
     805      199620 :       end subroutine update_snow_bgc
     806             : 
     807             : !=======================================================================
     808             : 
     809             :       end module icepack_aerosol
     810             : 
     811             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd