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

Generated by: LCOV version 1.14-6-g40580cd