Skip to content

Commit

Permalink
ICEPACK update
Browse files Browse the repository at this point in the history 
Update ICEPACK 1.2.4 to 1.3.4
but it needs further work, and other modules need to be updated as well
  • Loading branch information
@wangqian2284
wangqian2284 committed Sep 26, 2023
1 parent ddefd50 commit 4ba9298
Show file tree
Hide file tree
Showing 31 changed files with 6,184 additions and 4,936 deletions.
86 changes: 41 additions & 45 deletions src/Multi_ice/icepack_aerosol.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ module icepack_aerosol
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted

use icepack_zbgc_shared, only: kscavz
use icepack_zbgc_shared, only: kscavz

implicit none

Expand Down Expand Up @@ -59,14 +59,12 @@ subroutine update_aerosol(dt, &
aicen, & ! ice area fraction
aice_old, & ! values prior to thermodynamic changes
vice_old, &
vsno_old
vsno_old

real (kind=dbl_kind), dimension(:), &
intent(in) :: &
real (kind=dbl_kind), dimension(:), intent(in) :: &
faero_atm ! aerosol deposition rate (W/m^2 s)

real (kind=dbl_kind), dimension(:), &
intent(inout) :: &
real (kind=dbl_kind), dimension(:), intent(inout) :: &
faero_ocn ! aerosol flux to ocean (W/m^2 s)

real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
Expand Down Expand Up @@ -101,27 +99,27 @@ subroutine update_aerosol(dt, &

! echmod: this assumes max_aero=6
data kscav / .03_dbl_kind, .20_dbl_kind, .02_dbl_kind, &
.02_dbl_kind, .01_dbl_kind, .01_dbl_kind /
.02_dbl_kind, .01_dbl_kind, .01_dbl_kind /
data kscavsi / .03_dbl_kind, .20_dbl_kind, .02_dbl_kind, &
.02_dbl_kind, .01_dbl_kind, .01_dbl_kind /
.02_dbl_kind, .01_dbl_kind, .01_dbl_kind /

character(len=*),parameter :: subname='(update_aerosol)'

!-------------------------------------------------------------------
! initialize
!-------------------------------------------------------------------
focn_old(:) = faero_ocn(:)

hs_old = vsno_old/aice_old
hi_old = vice_old/aice_old
hslyr_old = hs_old/real(nslyr,kind=dbl_kind)
hilyr_old = hi_old/real(nilyr,kind=dbl_kind)

dzssl = min(hslyr_old/c2, hs_ssl)
dzssli = min(hilyr_old/c2, hi_ssl)
dzint = hs_old - dzssl
dzinti = hi_old - dzssli

if (aicen > c0) then
ar = c1/aicen
else ! ice disappeared during time step
Expand All @@ -146,7 +144,7 @@ subroutine update_aerosol(dt, &
aerotot0(k) = aerosno(k,2) + aerosno(k,1) &
+ aeroice(k,2) + aeroice(k,1)
enddo

!-------------------------------------------------------------------
! evaporation
!-------------------------------------------------------------------
Expand Down Expand Up @@ -212,7 +210,7 @@ subroutine update_aerosol(dt, &
aeroice(k,2) = aeroice(k,2) - sloss2
faero_ocn(k) = faero_ocn(k) + (sloss1+sloss2)/dt
enddo

dzinti = dzinti + min(dzssli+dhi_meltt, c0)
dzssli = max(dzssli+dhi_meltt, c0)
if (dzssli <= puny) then ! ssl ice melts away
Expand Down Expand Up @@ -251,7 +249,7 @@ subroutine update_aerosol(dt, &
enddo

dzssli = dzssli + min(dzinti+dhi_meltb, c0)
dzinti = max(dzinti+dhi_meltb, c0)
dzinti = max(dzinti+dhi_meltb, c0)
endif

!-------------------------------------------------------------------
Expand Down Expand Up @@ -292,7 +290,7 @@ subroutine update_aerosol(dt, &
else
hs = c0
endif
if (hs > hs_min) then ! should this really be hs_min or 0?
if (hs > hs_min) then ! should this really be hs_min or 0?
! should use same hs_min value as in radiation
do k=1,n_aero
aerosno(k,1) = aerosno(k,1) &
Expand Down Expand Up @@ -333,15 +331,15 @@ subroutine update_aerosol(dt, &
aeroice(k,1) = c0
enddo
endif

if (dzinti <= puny) then ! nothing in Ice INT
do k = 1, n_aero
faero_ocn(k) = faero_ocn(k) &
+ (aeroice(k,1)+aeroice(k,2))/dt
aeroice(k,:)=c0
enddo
endif

hslyr = hs/real(nslyr,kind=dbl_kind)
hilyr = hi/real(nilyr,kind=dbl_kind)
dzssl_new = min(hslyr/c2, hs_ssl)
Expand All @@ -356,32 +354,32 @@ subroutine update_aerosol(dt, &
dznew = max(dzssl_new-dzssl, c0)
if (dzint > puny) &
sloss1 = sloss1 + aerosno(k,2)*dznew/dzint
aerosno(k,1) = aerosno(k,1) + sloss1
aerosno(k,1) = aerosno(k,1) + sloss1
aerosno(k,2) = aerosno(k,2) - sloss1
enddo
else
aeroice(:,1) = aeroice(:,1) &
+ aerosno(:,1) + aerosno(:,2)
aerosno(:,:) = c0
endif

if (vicen > puny) then ! may want a limit on hi instead?
do k = 1, n_aero
sloss2 = c0
dznew = min(dzssli_new-dzssli, c0)
if (dzssli > puny) &
if (dzssli > puny) &
sloss2 = dznew*aeroice(k,1)/dzssli
dznew = max(dzssli_new-dzssli, c0)
if (dzinti > puny) &
if (dzinti > puny) &
sloss2 = sloss2 + aeroice(k,2)*dznew/dzinti
aeroice(k,1) = aeroice(k,1) + sloss2
aeroice(k,1) = aeroice(k,1) + sloss2
aeroice(k,2) = aeroice(k,2) - sloss2
enddo
else
faero_ocn(:) = faero_ocn(:) + (aeroice(:,1)+aeroice(:,2))/dt
aeroice(:,:) = c0
endif

!-------------------------------------------------------------------
! check conservation
!-------------------------------------------------------------------
Expand All @@ -406,7 +404,7 @@ subroutine update_aerosol(dt, &
! check for negative values
!-------------------------------------------------------------------

!echmod: note that this does not test or fix all aero tracers
!echmod: note that this does not test or fix all aero tracers
if (aeroice(1,1) < -puny .or. &
aeroice(1,2) < -puny .or. &
aerosno(1,1) < -puny .or. &
Expand Down Expand Up @@ -449,7 +447,7 @@ subroutine update_snow_bgc (dt, nblyr, &
ntrcr ! number of tracers

integer (kind=int_kind), dimension (nbtrcr), intent(in) :: &
bio_index
bio_index

real (kind=dbl_kind), intent(in) :: &
dt ! time step
Expand All @@ -466,19 +464,17 @@ subroutine update_snow_bgc (dt, nblyr, &
aicen, & ! ice area fraction
aice_old, & ! values prior to thermodynamic changes
vice_old, &
vsno_old
vsno_old

real (kind=dbl_kind),dimension(nbtrcr), intent(inout) :: &
zbgc_snow, & ! aerosol contribution from snow to ice
zbgc_atm, & ! and atm to ice concentration * volume (kg or mmol/m^3*m)
flux_bio ! total ocean tracer flux (mmol/m^2/s)

real (kind=dbl_kind), dimension(nbtrcr), &
intent(in) :: &
real (kind=dbl_kind), dimension(nbtrcr), intent(in) :: &
flux_bio_atm ! aerosol deposition rate (kg or mmol/m^2 s)

real (kind=dbl_kind), dimension(ntrcr), &
intent(inout) :: &
real (kind=dbl_kind), dimension(ntrcr), intent(inout) :: &
trcrn ! ice/snow tracer array

! local variables
Expand Down Expand Up @@ -547,21 +543,21 @@ subroutine update_snow_bgc (dt, nblyr, &
do k=1,nbtrcr
flux_bio(k) = flux_bio(k) + &
(trcrn(bio_index(k)+ nblyr+1)*dzssl+ &
trcrn(bio_index(k)+ nblyr+2)*dzint)/dt
trcrn(bio_index(k)+ nblyr+2)*dzint)/dt
trcrn(bio_index(k) + nblyr+1) = c0
trcrn(bio_index(k) + nblyr+2) = c0
zbgc_atm(k) = zbgc_atm(k) &
+ flux_bio_atm(k)*dt
+ flux_bio_atm(k)*dt
enddo

else
else

do k=1,nbtrcr
flux_bio_o(k) = flux_bio(k)
aerosno (k,1) = trcrn(bio_index(k)+ nblyr+1) * dzssl
aerosno (k,2) = trcrn(bio_index(k)+ nblyr+2) * dzint
aerosno0(k,:) = aerosno(k,:)
aerotot0(k) = aerosno(k,2) + aerosno(k,1)
aerotot0(k) = aerosno(k,2) + aerosno(k,1)
enddo

!-------------------------------------------------------------------
Expand All @@ -586,7 +582,7 @@ subroutine update_snow_bgc (dt, nblyr, &
*max(-dhs_melts-dzssl,c0)/dzint
aerosno(k,2) = aerosno(k,2) - sloss2
zbgc_snow(k) = zbgc_snow(k) + (sloss1+sloss2)
enddo !
enddo !

! update snow thickness
dzint=dzint+min(dzssl+dhs_melts, c0)
Expand All @@ -608,7 +604,7 @@ subroutine update_snow_bgc (dt, nblyr, &
!-------------------------------------------------------------------
! snowfall
!-------------------------------------------------------------------
if (fsnow > c0) dzssl = dzssl + fsnow/rhos*dt
if (fsnow > c0) dzssl = dzssl + fsnow/rhos*dt

!-------------------------------------------------------------------
! snow-ice formation
Expand Down Expand Up @@ -640,13 +636,13 @@ subroutine update_snow_bgc (dt, nblyr, &
else
hs = c0
endif
if (hs >= hs_min) then !should this really be hs_min or 0?
if (hs >= hs_min) then !should this really be hs_min or 0?
! should use same hs_min value as in radiation
do k=1,nbtrcr
aerosno(k,1) = aerosno(k,1) &
+ flux_bio_atm(k)*dt
enddo
else
else
do k=1,nbtrcr
zbgc_atm(k) = zbgc_atm(k) &
+ flux_bio_atm(k)*dt
Expand Down Expand Up @@ -678,7 +674,7 @@ subroutine update_snow_bgc (dt, nblyr, &
dzssl_new = min(hslyr/c2, hs_ssl)
dzint_new = hs - dzssl_new

if (hs > hs_min) then !should this really be hs_min or 0?
if (hs > hs_min) then !should this really be hs_min or 0?
do k = 1, nbtrcr
dznew = min(dzssl_new-dzssl, c0)
sloss1 = c0
Expand All @@ -687,7 +683,7 @@ subroutine update_snow_bgc (dt, nblyr, &
dznew = max(dzssl_new-dzssl, c0)
if (dzint > puny) &
sloss1 = sloss1 + aerosno(k,2)*dznew/dzint
aerosno(k,1) = aerosno(k,1) + sloss1
aerosno(k,1) = aerosno(k,1) + sloss1
aerosno(k,2) = aerosno(k,2) - sloss1
enddo
else
Expand All @@ -701,11 +697,11 @@ subroutine update_snow_bgc (dt, nblyr, &
!-------------------------------------------------------------------
do k = 1, nbtrcr
aerotot(k) = aerosno(k,2) + aerosno(k,1) &
+ zbgc_snow(k) + zbgc_atm(k)
+ zbgc_snow(k) + zbgc_atm(k)
aero_cons(k) = aerotot(k)-aerotot0(k) &
- ( flux_bio_atm(k) &
- ( flux_bio_atm(k) &
- (flux_bio(k)-flux_bio_o(k))) * dt
if (aero_cons(k) > puny .or. zbgc_snow(k) + zbgc_atm(k) < c0) then
if (aero_cons(k) > puny .or. zbgc_snow(k) + zbgc_atm(k) < c0) then
write(warnstr,*) subname, 'Conservation failure: aerosols in snow'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, 'test aerosol 1'
Expand Down Expand Up @@ -735,7 +731,7 @@ subroutine update_snow_bgc (dt, nblyr, &
!-------------------------------------------------------------------
if (vsnon > puny) then
do k = 1,nbtrcr
trcrn(bio_index(k)+nblyr+1)=aerosno(k,1)/dzssl_new
trcrn(bio_index(k)+nblyr+1)=aerosno(k,1)/dzssl_new
trcrn(bio_index(k)+nblyr+2)=aerosno(k,2)/dzint_new
enddo
else
Expand Down
5 changes: 2 additions & 3 deletions src/Multi_ice/icepack_age.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,13 +26,12 @@ subroutine increment_age (dt, iage)
real (kind=dbl_kind), intent(in) :: &
dt ! time step

real (kind=dbl_kind), &
intent(inout) :: &
real (kind=dbl_kind), intent(inout) :: &
iage

character(len=*),parameter :: subname='(increment_age)'

iage = iage + dt
iage = iage + dt

end subroutine increment_age

Expand Down
Loading

0 comments on commit 4ba9298

Please sign in to comment.