!=======================================================================
!
!BOP
!
! !MODULE: ice_aerosol - Aerosol tracer within sea ice
!
! !DESCRIPTION:
!
! !REVISION HISTORY:
! SVN:$$
!
! authors Marika Holland, NCAR
! David Bailey, NCAR
!
! !INTERFACE:
!
module ice_aerosol 4,6
!
! !USES:
!
use ice_kinds_mod
use ice_constants
use ice_fileunits
use ice_restart, only: lenstr, restart_dir, restart_file, &
pointer_file, runtype
use ice_communicate, only: my_task, master_task
use ice_exit, only: abort_ice
!
!EOP
!
implicit none
logical (kind=log_kind) :: &
restart_aero ! if .true., read aerosol tracer restart file
!=======================================================================
contains
!=======================================================================
!BOP
!
! !ROUTINE: init_aerosol
!
! !DESCRIPTION:
!
! Initialize ice aerosol tracer (call prior to reading restart data)
!
! !REVISION HISTORY: same as module
!
! !INTERFACE:
!
subroutine init_aerosol 1,3
!
! !USES:
!
use ice_state, only: filename_aero
!
!EOP
!
if (trim(filename_aero) /= 'none') restart_aero = .true.
if (restart_aero) then
if (trim(runtype) == 'continue') then
call read_restart_aero
else
call read_restart_aero(filename_aero)
endif
endif
end subroutine init_aerosol
!=======================================================================
!BOP
!
! !ROUTINE: update_aerosol
!
! !DESCRIPTION:
!
! Increase aerosol in ice or snow surface due to deposition
!
! !REVISION HISTORY: same as module
!
! !INTERFACE:
!
subroutine update_aerosol (nx_block, ny_block, & 1,2
dt, icells, &
indxi, indxj, &
meltt, melts, &
meltb, congel, &
snoice, &
fsnow, &
trcrn, &
aice_old, &
vice_old, vsno_old, &
vicen, vsnon, aicen, &
faero, fsoot)
!
! !USES:
!
use ice_domain_size, only: max_ntrcr, nilyr, nslyr, n_aero, n_aeromx
use ice_state, only: nt_aero
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
icells ! number of cells with ice present
integer (kind=int_kind), dimension (nx_block*ny_block), &
intent(in) :: &
indxi, indxj ! compressed indices for cells with ice
real (kind=dbl_kind), intent(in) :: &
dt ! time step
real (kind=dbl_kind), dimension(nx_block,ny_block), &
intent(in) :: &
meltt, &
melts, &
meltb, &
congel, &
snoice, &
fsnow, &
vicen, &
vsnon, &
aicen, &
aice_old, &
vice_old, &
vsno_old
real (kind=dbl_kind), dimension(nx_block,ny_block,n_aeromx), &
intent(in) :: &
faero
real (kind=dbl_kind), dimension(nx_block,ny_block,n_aeromx), &
intent(inout) :: &
fsoot
real (kind=dbl_kind), dimension(nx_block,ny_block,max_ntrcr), &
intent(inout) :: &
trcrn
!
! local variables
!
integer (kind=int_kind) :: i, j, ij, k
integer (kind=int_kind) :: n ! print_points
!
real (kind=dbl_kind), dimension(icells) :: &
dzssl, &
dzint, &
dzssli, &
dzinti
real (kind=dbl_kind), dimension(icells) :: &
dhs_evap, dhi_evap, &
dhs_melts, dhs_snoice, dhi_meltt, dhi_snoice, &
dhi_congel, dhi_meltb
real (kind=dbl_kind), dimension(icells,n_aeromx) :: &
aerotot, aerotot0 ! for diagnostics
real (kind=dbl_kind) :: &
dzssl_new, &
dzint_new, &
dzssli_new, &
dzinti_new, &
dznew
real (kind=dbl_kind), dimension(nx_block,ny_block,n_aeromx,2) :: &
aerosno, aeroice, &
aerosno0, aeroice0 ! for diagnostic prints
real (kind=dbl_kind), dimension(nx_block,ny_block,n_aeromx) :: &
fsoot_old
real (kind=dbl_kind) :: &
hs_old, hi_old, hslyr_old, hilyr_old, dhs, dhi, hs, hi, &
hslyr, hilyr, sloss1, sloss2
real (kind=dbl_kind), dimension(n_aeromx) :: &
kscav, kscavsi
!MH These need to be the same as in the DE code. Put in a common place?
real (kind=dbl_kind) :: &
hi_ssl, hs_ssl
data hs_ssl / .040_dbl_kind /
data hi_ssl / .050_dbl_kind /
data kscav / .03_dbl_kind, .20_dbl_kind,&
.02_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 /
aerosno(:,:,:,:) = c0
aeroice(:,:,:,:) = c0
aerosno0(:,:,:,:) = c0
aeroice0(:,:,:,:) = c0
fsoot_old(:,:,:) = fsoot(:,:,:)
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
hs_old=vsno_old(i,j)/aice_old(i,j)
hi_old=vice_old(i,j)/aice_old(i,j)
hslyr_old=hs_old/real(nslyr,kind=dbl_kind)
hilyr_old=hi_old/real(nilyr,kind=dbl_kind)
dzssl(ij)=min(hslyr_old/c2,hs_ssl)
dzint(ij)=hs_old-dzssl(ij)
dzssli(ij)=min(hilyr_old/c2,hi_ssl)
dzinti(ij)=hi_old-dzssli(ij)
if (aicen(i,j) > c0) then
hs = vsnon(i,j)/aicen(i,j)
hi = vicen(i,j)/aicen(i,j)
dhs_melts(ij)=-melts(i,j)/aicen(i,j)
dhi_snoice(ij)=snoice(i,j)/aicen(i,j)
dhs_snoice(ij)=dhi_snoice(ij)*rhoi/rhos
dhi_meltt(ij)=-meltt(i,j)/aicen(i,j)
dhi_meltb(ij)=-meltb(i,j)/aicen(i,j)
dhi_congel(ij)=congel(i,j)/aicen(i,j)
else
hs = vsnon(i,j)/aice_old(i,j)
hi = vicen(i,j)/aice_old(i,j)
dhs_melts(ij)=-melts(i,j)/aice_old(i,j)
dhi_snoice(ij)=snoice(i,j)/aice_old(i,j)
dhs_snoice(ij)=dhi_snoice(ij)*rhoi/rhos
dhi_meltt(ij)=-meltt(i,j)/aice_old(i,j)
dhi_meltb(ij)=-meltb(i,j)/aice_old(i,j)
dhi_congel(ij)=congel(i,j)/aice_old(i,j)
endif
dhs_evap(ij)=hs-(hs_old+dhs_melts(ij)-dhs_snoice(ij)+&
fsnow(i,j)/rhos*dt)
dhi_evap(ij)=hi-(hi_old+dhi_meltt(ij)+dhi_meltb(ij)+ &
dhi_congel(ij)+dhi_snoice(ij))
enddo
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
do k=1,n_aero
aerosno(i,j,k,:)=&
trcrn(i,j,nt_aero+(k-1)*4 :nt_aero+(k-1)*4+1)*vsno_old(i,j) ! aerosol in snow
aeroice(i,j,k,:)=&
trcrn(i,j,nt_aero+(k-1)*4+2:nt_aero+(k-1)*4+3)*vice_old(i,j) ! aerosol in ice
aerosno0(i,j,k,:)=aerosno(i,j,k,:)
aeroice0(i,j,k,:)=aeroice(i,j,k,:)
aerotot0(ij,k)=aerosno(i,j,k,2)+aerosno(i,j,k,1) &
+aeroice(i,j,k,2)+aeroice(i,j,k,1)
enddo
enddo
! apply evaporation
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
dzint(ij)=dzint(ij) + min(dzssl(ij)+dhs_evap(ij),c0)
dzssl(ij)=max(dzssl(ij)+dhs_evap(ij),c0)
dzinti(ij)=dzinti(ij) + min(dzssli(ij)+dhi_evap(ij),c0)
dzssli(ij)=max(dzssli(ij)+dhi_evap(ij),c0)
enddo
! basal ice growth
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
dzinti(ij)=dzinti(ij)+dhi_congel(ij)
enddo
! surface snow melt
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (-dhs_melts(ij) > puny) then
do k=1,n_aero
sloss1=c0
sloss2=c0
if (dzssl(ij) > puny) &
sloss1=kscav(k)*aerosno(i,j,k,1) &
*min(-dhs_melts(ij),dzssl(ij))/dzssl(ij)
aerosno(i,j,k,1)=aerosno(i,j,k,1)-sloss1
if (dzint(ij) > puny) &
sloss2=kscav(k)*aerosno(i,j,k,2) &
*max(-dhs_melts(ij)-dzssl(ij),c0)/dzint(ij)
aerosno(i,j,k,2)=aerosno(i,j,k,2)-sloss2
fsoot(i,j,k)=fsoot(i,j,k)+(sloss1+sloss2)/dt
enddo ! n_aero
! update snow thickness
dzint(ij)=dzint(ij)+min(dzssl(ij)+dhs_melts(ij),c0)
dzssl(ij)=max(dzssl(ij)+dhs_melts(ij),c0)
if ( dzssl(ij) <= puny ) then ! ssl melts away
aerosno(i,j,:,2)=aerosno(i,j,:,1)+aerosno(i,j,:,2)
aerosno(i,j,:,1)=c0
dzssl(ij)=max(dzssl(ij),c0)
endif
if (dzint(ij) <= puny ) then ! all snow melts away
aeroice(i,j,:,1)=&
aeroice(i,j,:,1)+aerosno(i,j,:,1)+aerosno(i,j,:,2)
aerosno(i,j,:,:)=c0
dzint(ij)=max(dzint(ij),c0)
endif
endif
enddo
! surface ice melt
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (-dhi_meltt(ij) > puny) then
do k=1,n_aero
sloss1=c0
sloss2=c0
if (dzssli(ij) > puny) &
sloss1=kscav(k)*aeroice(i,j,k,1) &
*min(-dhi_meltt(ij),dzssli(ij))/dzssli(ij)
aeroice(i,j,k,1)=aeroice(i,j,k,1)-sloss1
if (dzinti(ij) > puny) &
sloss2=kscav(k)*aeroice(i,j,k,2) &
*max(-dhi_meltt(ij)-dzssli(ij),c0)/dzinti(ij)
aeroice(i,j,k,2)=aeroice(i,j,k,2)-sloss2
fsoot(i,j,k)=fsoot(i,j,k)+(sloss1+sloss2)/dt
enddo
dzinti(ij)=dzinti(ij)+min(dzssli(ij)+dhi_meltt(ij),c0)
dzssli(ij)=max(dzssli(ij)+dhi_meltt(ij),c0)
if (dzssli(ij) <= puny) then ! ssl ice melts away
do k=1,n_aero
aeroice(i,j,k,2)=aeroice(i,j,k,1)+aeroice(i,j,k,2)
aeroice(i,j,k,1)=c0
enddo
dzssli(ij)=max(dzssli(ij),c0)
endif
if (dzinti(ij) <= puny) then ! all ice melts away
do k=1,n_aero
fsoot(i,j,k)=fsoot(i,j,k) &
+(aeroice(i,j,k,1)+aeroice(i,j,k,2))/dt
aeroice(i,j,k,:)=c0
enddo
dzinti(ij)=max(dzinti(ij),c0)
endif
endif
enddo
! basal ice melt. Assume all soot lost in basal melt
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (-dhi_meltb(ij) > puny) then
do k=1,n_aero
sloss1=c0
sloss2=c0
if (dzssli(ij) > puny) &
sloss1=max(-dhi_meltb(ij)-dzinti(ij),c0) &
*aeroice(i,j,k,1)/dzssli(ij)
aeroice(i,j,k,1)=aeroice(i,j,k,1)-sloss1
if (dzinti(ij) > puny) &
sloss2=min(-dhi_meltb(ij),dzinti(ij)) &
*aeroice(i,j,k,2)/dzinti(ij)
aeroice(i,j,k,2)=aeroice(i,j,k,2)-sloss2
fsoot(i,j,k)=fsoot(i,j,k)+(sloss1+sloss2)/dt
enddo
dzssli(ij) = dzssli(ij)+min(dzinti(ij)+dhi_meltb(ij), c0)
dzinti(ij) = max(dzinti(ij)+dhi_meltb(ij), c0)
endif
enddo
! snowfall
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (fsnow(i,j) > c0) &
dzssl(ij)=dzssl(ij)+fsnow(i,j)/rhos*dt
enddo
! snoice formation
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (dhs_snoice(ij) > puny) then
do k=1,n_aero
sloss1=c0
sloss2=c0
if (dzint(ij) > puny) &
sloss2 = min(dhs_snoice(ij),dzint(ij)) &
*aerosno(i,j,k,2)/dzint(ij)
aerosno(i,j,k,2) = aerosno(i,j,k,2) - sloss2
if (dzssl(ij) > puny) &
sloss1 = max(dhs_snoice(ij)-dzint(ij),c0) &
*aerosno(i,j,k,1)/dzssl(ij)
aerosno(i,j,k,1) = aerosno(i,j,k,1) - sloss1
aeroice(i,j,k,1) = aeroice(i,j,k,1) &
+ (c1-kscavsi(k))*(sloss2+sloss1)
fsoot(i,j,k)=fsoot(i,j,k)+kscavsi(k)*(sloss2+sloss1)/dt
enddo
dzssl(ij)=dzssl(ij)-max(dhs_snoice(ij)-dzint(ij),c0)
dzint(ij)=max(dzint(ij)-dhs_snoice(ij),c0)
dzssli(ij)=dzssli(ij)+dhi_snoice(ij)
endif
enddo
! aerosol deposition
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (aicen(i,j) > c0) then
hs = vsnon(i,j) / aicen(i,j)
else
hs = c0
endif
if (hs > hsmin) then ! should this really be hsmin or 0?
! should use same hsmin value as in radiation
do k=1,n_aero
aerosno(i,j,k,1)=aerosno(i,j,k,1) &
+ faero(i,j,k)*dt*aicen(i,j)
enddo
else
do k=1,n_aero
aeroice(i,j,k,1)=aeroice(i,j,k,1) &
+ faero(i,j,k)*dt*aicen(i,j)
enddo
endif
enddo
! redistribute aerosol within vertical layers
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (aicen(i,j) > c0) then
hs = vsnon(i,j) / aicen(i,j) ! new snow thickness
hi = vicen(i,j) / aicen(i,j) ! new ice thickness
else
hs = c0
hi = c0
endif
if (dzssl(ij) <= puny) then ! nothing in SSL
do k=1,n_aero
aerosno(i,j,k,2)=aerosno(i,j,k,2)+aerosno(i,j,k,1)
aerosno(i,j,k,1)=c0
enddo
endif
if (dzint(ij) <= puny) then ! nothing in Snow Int
do k=1,n_aero
aeroice(i,j,k,1)=aeroice(i,j,k,1)+aerosno(i,j,k,2)
aerosno(i,j,k,2)=c0
enddo
endif
if (dzssli(ij) <= puny) then ! nothing in Ice SSL
do k=1,n_aero
aeroice(i,j,k,2)=aeroice(i,j,k,2)+aeroice(i,j,k,1)
aeroice(i,j,k,1)=c0
enddo
endif
if (dzinti(ij) <= puny) then ! nothing in Ice INT
do k=1,n_aero
fsoot(i,j,k)=fsoot(i,j,k)+&
(aeroice(i,j,k,1)+aeroice(i,j,k,2))/dt
aeroice(i,j,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) ! ssl for snow
dzint_new=hs-dzssl_new
dzssli_new=min(hilyr/c2,hi_ssl) ! ssl for ice
dzinti_new=hi-dzssli_new
if (hs > hsmin) then
do k=1,n_aero
dznew=min(dzssl_new-dzssl(ij),c0)
sloss1=c0
if (dzssl(ij) > puny) &
sloss1=dznew*aerosno(i,j,k,1)/dzssl(ij) ! not neccesarily a loss term
dznew=max(dzssl_new-dzssl(ij),c0)
if (dzint(ij) > puny) &
sloss1=sloss1+aerosno(i,j,k,2)*dznew/dzint(ij) ! not really a loss term
aerosno(i,j,k,1) =aerosno(i,j,k,1)+sloss1
aerosno(i,j,k,2) =aerosno(i,j,k,2)-sloss1
enddo
else
aeroice(i,j,:,1)=aeroice(i,j,:,1) &
+aerosno(i,j,:,1)+aerosno(i,j,:,2)
aerosno(i,j,:,:) = c0
endif
if (vicen(i,j) > puny) then ! may want a limit on hi instead?
do k=1,n_aero
sloss2=c0
dznew=min(dzssli_new-dzssli(ij),c0)
if (dzssli(ij) > puny) &
sloss2=dznew*aeroice(i,j,k,1)/dzssli(ij)
dznew=max(dzssli_new-dzssli(ij),c0)
if (dzinti(ij) > puny) &
sloss2=sloss2+aeroice(i,j,k,2)*dznew/dzinti(ij) ! not really a loss term
aeroice(i,j,k,1) =aeroice(i,j,k,1)+sloss2
aeroice(i,j,k,2) =aeroice(i,j,k,2)-sloss2
enddo
else
fsoot(i,j,:)=fsoot(i,j,:)+(aeroice(i,j,:,1)+aeroice(i,j,:,2))/dt
aeroice(i,j,:,:) = c0
endif
do k=1,n_aero
aerotot(ij,k)=aerosno(i,j,k,2)+aerosno(i,j,k,1) &
+aeroice(i,j,k,2)+aeroice(i,j,k,1)
if ( ( (aerotot(ij,k)-aerotot0(ij,k)) &
- ( faero(i,j,k)*aicen(i,j) &
- (fsoot(i,j,k)-fsoot_old(i,j,k)) )*dt ) > 0.00001) then
write(nu_diag,*) 'aerosol tracer: ',k
write(nu_diag,*) 'aerotot-aerotot0 ',aerotot(ij,k)-aerotot0(ij,k)
write(nu_diag,*) 'faero-fsoot ',faero(i,j,k)*aicen(i,j)*dt &
-(fsoot(i,j,k)-fsoot_old(i,j,k))*dt
endif
enddo
enddo
! reload tracers
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (vicen(i,j) > puny) &
aeroice(i,j,:,:)=aeroice(i,j,:,:)/vicen(i,j)
if (vsnon(i,j) > puny) &
aerosno(i,j,:,:)=aerosno(i,j,:,:)/vsnon(i,j)
do k=1,n_aero
do n=1,2
trcrn(i,j,nt_aero+(k-1)*4+n-1)=aerosno(i,j,k,n)
trcrn(i,j,nt_aero+(k-1)*4+n+1)=aeroice(i,j,k,n)
enddo
! do n=1,4
! if (trcrn(i,j,nt_aero+(k-1)*4+n-1) < puny) then
! fsoot(i,j,k)=fsoot(i,j,k)+ &
! trcrn(i,j,nt_aero+(k-1)*4+n-1)/dt
! trcrn(i,j,nt_aero+(k-1)*4+n-1)=c0
! endif
! enddo
enddo
enddo
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (trcrn(i,j,nt_aero) < -puny .or. trcrn(i,j,nt_aero+1) < -puny &
.or. trcrn(i,j,nt_aero+2) < -puny .or. trcrn(i,j,nt_aero+3) < -puny) then
write(nu_diag,*) 'MH aerosol negative in aerosol code'
write(nu_diag,*) 'MH INT neg in aerosol my_task = ',&
my_task &
,' printing point = ',n &
,' i and j = ',i,j
write(nu_diag,*) 'MH Int Neg aero snowssl= ',aerosno0(i,j,1,1)
write(nu_diag,*) 'MH Int Neg aero new snowssl= ',aerosno(i,j,1,1)
write(nu_diag,*) 'MH Int Neg aero snowint= ',aerosno0(i,j,1,2)
write(nu_diag,*) 'MH Int Neg aero new snowint= ',aerosno(i,j,1,2)
write(nu_diag,*) 'MH Int Neg aero ice_ssl= ',aeroice0(i,j,1,1)
write(nu_diag,*) 'MH Int Neg aero new ice_ssl= ',aeroice(i,j,1,1)
write(nu_diag,*) 'MH Int Neg aero ice_int= ',aeroice0(i,j,1,2)
write(nu_diag,*) 'MH Int Neg aero new ice_int= ',aeroice(i,j,1,2)
write(nu_diag,*) 'MH Int Neg aero aicen= ',aicen(i,j)
write(nu_diag,*) 'MH Int Neg aero vicen= ',vicen(i,j)
write(nu_diag,*) 'MH Int Neg aero vsnon= ',vsnon(i,j)
write(nu_diag,*) 'MH Int Neg aero viceold= ',vice_old(i,j)
write(nu_diag,*) 'MH Int Neg aero vsnoold= ',vsno_old(i,j)
write(nu_diag,*) 'MH Int Neg aero melts= ',melts(i,j)
write(nu_diag,*) 'MH Int Neg aero meltt= ',meltt(i,j)
write(nu_diag,*) 'MH Int Neg aero meltb= ',meltb(i,j)
write(nu_diag,*) 'MH Int Neg aero congel= ',congel(i,j)
write(nu_diag,*) 'MH Int Neg aero snoice= ',snoice(i,j)
write(nu_diag,*) 'MH Int Neg aero evap sno?= ',dhs_evap(ij)
write(nu_diag,*) 'MH Int Neg aero evap ice?= ',dhi_evap(ij)
write(nu_diag,*) 'MH Int Neg aero fsnow= ',fsnow(i,j)
write(nu_diag,*) 'MH Int Neg aero faero= ',faero(i,j,1)
write(nu_diag,*) 'MH Int Neg aero fsoot= ',fsoot(i,j,1)
trcrn(i,j,nt_aero)=max(trcrn(i,j,nt_aero),c0)
trcrn(i,j,nt_aero+1)=max(trcrn(i,j,nt_aero+1),c0)
trcrn(i,j,nt_aero+2)=max(trcrn(i,j,nt_aero+2),c0)
trcrn(i,j,nt_aero+3)=max(trcrn(i,j,nt_aero+3),c0)
endif
enddo
end subroutine update_aerosol
!=======================================================================
!---! these subroutines write/read Fortran unformatted data files ..
!=======================================================================
!
!BOP
!
! !IROUTINE: write_restart_aero - dumps all fields required for restart
!
! !INTERFACE:
!
subroutine write_restart_aero(filename_spec) 1,10
!
! !DESCRIPTION:
!
! Dumps all values needed for restarting
!
! !REVISION HISTORY:
!
! authors Elizabeth Hunke, LANL (original version)
! David Bailey, NCAR
! Marika Holland, NCAR
!
! !USES:
!
use ice_domain_size
use ice_calendar, only: sec, month, mday, nyr, istep1, &
time, time_forc, idate, year_init
use ice_state
use ice_read_write
use ice_restart, only: lenstr, restart_dir, restart_file, pointer_file
!
! !INPUT/OUTPUT PARAMETERS:
!
character(len=char_len_long), intent(in), optional :: filename_spec
!EOP
!
integer (kind=int_kind) :: &
i, j, k, n, it, iblk, & ! counting indices
iyear, imonth, iday ! year, month, day
character(len=char_len_long) :: filename
logical (kind=log_kind) :: diag
! construct path/file
if (present(filename_spec)) then
filename = trim(filename_spec)
else
iyear = nyr + year_init - 1
imonth = month
iday = mday
write(filename,'(a,a,a,i4.4,a,i2.2,a,i2.2,a,i5.5)') &
restart_dir(1:lenstr(restart_dir)), &
restart_file(1:lenstr(restart_file)),'.aero.', &
iyear,'-',month,'-',mday,'-',sec
end if
! begin writing restart data
call ice_open(nu_dump_aero,filename,0)
if (my_task == master_task) then
write(nu_dump_aero) istep1,time,time_forc
write(nu_diag,*) 'Writing ',filename(1:lenstr(filename))
endif
diag = .true.
!-----------------------------------------------------------------
do k = 1, n_aero
do n = 1, ncat
call ice_write(nu_dump_aero,0,trcrn(:,:,nt_aero +(k-1)*4,n,:),'ruf8',diag)
call ice_write(nu_dump_aero,0,trcrn(:,:,nt_aero+1+(k-1)*4,n,:),'ruf8',diag)
call ice_write(nu_dump_aero,0,trcrn(:,:,nt_aero+2+(k-1)*4,n,:),'ruf8',diag)
call ice_write(nu_dump_aero,0,trcrn(:,:,nt_aero+3+(k-1)*4,n,:),'ruf8',diag)
enddo
enddo
if (my_task == master_task) close(nu_dump_aero)
end subroutine write_restart_aero
!=======================================================================
!BOP
!
! !IROUTINE: read_restart_aero - reads all fields required for restart
!
! !INTERFACE:
!
subroutine read_restart_aero(filename_spec) 2,12
!
! !DESCRIPTION:
!
! Reads all values needed for an ice aerosol restart
!
! !REVISION HISTORY:
!
! authors Elizabeth Hunke, LANL (original version)
! David Bailey, NCAR
! Marika Holland, NCAR
!
! !USES:
!
use ice_domain_size
use ice_calendar, only: sec, month, mday, nyr, istep1, &
time, time_forc, idate, year_init
use ice_state
use ice_read_write
use ice_restart, only: lenstr, restart_dir, restart_file, pointer_file
!
! !INPUT/OUTPUT PARAMETERS:
!
character(len=char_len_long), intent(in), optional :: filename_spec
!EOP
!
integer (kind=int_kind) :: &
i, j, k, n, it, iblk, & ! counting indices
iyear, imonth, iday ! year, month, day
character(len=char_len_long) :: &
filename, filename0, string1, string2
logical (kind=log_kind) :: &
diag
if (my_task == master_task) then
! reconstruct path/file
if (present(filename_spec)) then
filename = filename_spec
else
open(nu_rst_pointer,file=pointer_file)
read(nu_rst_pointer,'(a)') filename0
filename = trim(filename0)
close(nu_rst_pointer)
n = index(filename0,trim(restart_file))
if (n == 0) call abort_ice('soot restart: filename discrepancy')
string1 = trim(filename0(1:n-1))
string2 = trim(filename0(n+lenstr(restart_file):lenstr(filename0)))
write(filename,'(a,a,a,a)') &
string1(1:lenstr(string1)), &
restart_file(1:lenstr(restart_file)),'.aero', &
string2(1:lenstr(string2))
endif
endif ! master_task
call ice_open(nu_restart_aero,filename,0)
if (my_task == master_task) then
read(nu_restart_aero) istep1,time,time_forc
write(nu_diag,*) 'Reading ',filename(1:lenstr(filename))
endif
diag = .true.
!-----------------------------------------------------------------
do k = 1, n_aero
do n = 1, ncat
call ice_read(nu_restart_aero,0,trcrn(:,:,nt_aero +(k-1)*4,n,:),'ruf8',&
diag,field_type=field_type_scalar,field_loc=field_loc_center)
call ice_read(nu_restart_aero,0,trcrn(:,:,nt_aero+1+(k-1)*4,n,:),'ruf8',&
diag,field_type=field_type_scalar,field_loc=field_loc_center)
call ice_read(nu_restart_aero,0,trcrn(:,:,nt_aero+2+(k-1)*4,n,:),'ruf8',&
diag,field_type=field_type_scalar,field_loc=field_loc_center)
call ice_read(nu_restart_aero,0,trcrn(:,:,nt_aero+3+(k-1)*4,n,:),'ruf8',&
diag,field_type=field_type_scalar,field_loc=field_loc_center)
enddo
enddo
if (my_task == master_task) close(nu_restart_aero)
end subroutine read_restart_aero
!=======================================================================
end module ice_aerosol
!=======================================================================