-
Notifications
You must be signed in to change notification settings - Fork 0
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Dyninit #14
Dyninit #14
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -42,7 +42,8 @@ module ice_dyn_evp | |
use ice_constants, only: c0, p027, p055, p111, p166, & | ||
p222, p25, p333, p5, c1 | ||
use ice_dyn_shared, only: stepu, dyn_prep1, dyn_prep2, dyn_finish, & | ||
ndte, yield_curve, ecci, denom1, arlx1i, fcor_blk, uvel_init, vvel_init, & | ||
ndte, yield_curve, ecci, denom1, arlx1i, fcor_blk, fcorE_blk, fcorN_blk, & | ||
uvel_init, vvel_init, uvelE_init, vvelE_init, uvelN_init, vvelN_init, & | ||
seabed_stress_factor_LKD, seabed_stress_factor_prob, seabed_stress_method, & | ||
seabed_stress, Ktens, revp | ||
use ice_fileunits, only: nu_diag | ||
|
@@ -84,14 +85,22 @@ subroutine evp (dt) | |
strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, & | ||
strtltx, strtlty, strocnx, strocny, strintx, strinty, taubx, tauby, & | ||
strocnxT, strocnyT, strax, stray, & | ||
strairxN, strairyN, uocnN, vocnN, ss_tltxN, ss_tltyN, icenmask, fmN, & | ||
strtltxN, strtltyN, strocnxN, strocnyN, strintxN, strintyN, taubxN, taubyN, & | ||
straxN, strayN, TbN, & | ||
strairxE, strairyE, uocnE, vocnE, ss_tltxE, ss_tltyE, iceemask, fmE, & | ||
strtltxE, strtltyE, strocnxE, strocnyE, strintxE, strintyE, taubxE, taubyE, & | ||
straxE, strayE, TbE, & | ||
Tbu, hwater, & | ||
stressp_1, stressp_2, stressp_3, stressp_4, & | ||
stressm_1, stressm_2, stressm_3, stressm_4, & | ||
stress12_1, stress12_2, stress12_3, stress12_4 | ||
use ice_grid, only: tmask, umask, dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, & | ||
use ice_grid, only: tmask, umask, nmask, emask, dxt, dyt, & | ||
dxhy, dyhx, cxp, cyp, cxm, cym, & | ||
tarear, uarear, tinyarea, grid_average_X2Y, & | ||
grid_type, grid_system | ||
use ice_state, only: aice, vice, vsno, uvel, vvel, divu, shear, & | ||
use ice_state, only: aice, vice, vsno, uvel, vvel, uvelN, vvelN, & | ||
uvelE, vvelE, divu, shear, & | ||
aice_init, aice0, aicen, vicen, strength | ||
use ice_timers, only: timer_dynamics, timer_bound, & | ||
ice_timer_start, ice_timer_stop, timer_evp_1d, timer_evp_2d | ||
|
@@ -112,11 +121,17 @@ subroutine evp (dt) | |
|
||
integer (kind=int_kind), dimension(max_blocks) :: & | ||
icellt , & ! no. of cells where icetmask = 1 | ||
icelln , & ! no. of cells where icenmask = 1 | ||
icelle , & ! no. of cells where iceemask = 1 | ||
icellu ! no. of cells where iceumask = 1 | ||
|
||
integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks) :: & | ||
indxti , & ! compressed index in i-direction | ||
indxtj , & ! compressed index in j-direction | ||
indxei , & ! compressed index in i-direction | ||
indxej , & ! compressed index in j-direction | ||
indxni , & ! compressed index in i-direction | ||
indxnj , & ! compressed index in j-direction | ||
indxui , & ! compressed index in i-direction | ||
indxuj ! compressed index in j-direction | ||
|
||
|
@@ -130,6 +145,24 @@ subroutine evp (dt) | |
umass , & ! total mass of ice and snow (u grid) | ||
umassdti ! mass of U-cell/dte (kg/m^2 s) | ||
|
||
real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & | ||
waterxN , & ! for ocean stress calculation, x (m/s) | ||
wateryN , & ! for ocean stress calculation, y (m/s) | ||
forcexN , & ! work array: combined atm stress and ocn tilt, x | ||
forceyN , & ! work array: combined atm stress and ocn tilt, y | ||
aiN , & ! ice fraction on N-grid | ||
nmass , & ! total mass of ice and snow (N grid) | ||
nmassdti ! mass of N-cell/dte (kg/m^2 s) | ||
|
||
real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & | ||
waterxE , & ! for ocean stress calculation, x (m/s) | ||
wateryE , & ! for ocean stress calculation, y (m/s) | ||
forcexE , & ! work array: combined atm stress and ocn tilt, x | ||
forceyE , & ! work array: combined atm stress and ocn tilt, y | ||
aiE , & ! ice fraction on E-grid | ||
emass , & ! total mass of ice and snow (E grid) | ||
emassdti ! mass of E-cell/dte (kg/m^2 s) | ||
|
||
real (kind=dbl_kind), allocatable :: fld2(:,:,:,:) | ||
|
||
real (kind=dbl_kind), dimension(nx_block,ny_block,8):: & | ||
|
@@ -207,6 +240,12 @@ subroutine evp (dt) | |
strairx (:,:,iblk), strairy (:,:,iblk), & | ||
tmass (:,:,iblk), icetmask(:,:,iblk)) | ||
|
||
if (grid_system == 'CD') then | ||
strairxN(:,:,iblk) = strairxT(:,:,iblk) | ||
strairyN(:,:,iblk) = strairyT(:,:,iblk) | ||
strairxE(:,:,iblk) = strairxT(:,:,iblk) | ||
strairyE(:,:,iblk) = strairyT(:,:,iblk) | ||
endif | ||
enddo ! iblk | ||
!$OMP END PARALLEL DO | ||
|
||
|
@@ -222,6 +261,12 @@ subroutine evp (dt) | |
call grid_average_X2Y('T2UF',tmass,umass) | ||
call grid_average_X2Y('T2UF',aice_init, aiu) | ||
|
||
if (grid_system == 'CD') then | ||
call grid_average_X2Y('T2EF',tmass,emass) | ||
call grid_average_X2Y('T2EF',aice_init, aie) | ||
call grid_average_X2Y('T2NF',tmass,nmass) | ||
call grid_average_X2Y('T2NF',aice_init, aie) | ||
endif | ||
!---------------------------------------------------------------- | ||
! Set wind stress to values supplied via NEMO or other forcing | ||
! This wind stress is rotated on u grid and multiplied by aice | ||
|
@@ -243,6 +288,29 @@ subroutine evp (dt) | |
call grid_average_X2Y('T2UF',strairy) | ||
endif | ||
|
||
if (grid_system == 'CD') then | ||
|
||
if (.not. calc_strair) then | ||
strairxN(:,:,:) = strax(:,:,:) | ||
strairyN(:,:,:) = stray(:,:,:) | ||
strairxE(:,:,:) = strax(:,:,:) | ||
strairyE(:,:,:) = stray(:,:,:) | ||
else | ||
call ice_HaloUpdate (strairxN, halo_info, & | ||
field_loc_center, field_type_vector) | ||
call ice_HaloUpdate (strairyN, halo_info, & | ||
field_loc_center, field_type_vector) | ||
call ice_HaloUpdate (strairxE, halo_info, & | ||
field_loc_center, field_type_vector) | ||
call ice_HaloUpdate (strairyE, halo_info, & | ||
field_loc_center, field_type_vector) | ||
call grid_average_X2Y('T2NF',strairxN) | ||
call grid_average_X2Y('T2NF',strairyN) | ||
call grid_average_X2Y('T2EF',strairxE) | ||
call grid_average_X2Y('T2EF',strairyE) | ||
endif | ||
|
||
endif | ||
! tcraig, tcx, threading here leads to some non-reproducbile results and failures in icepack_ice_strength | ||
! need to do more debugging | ||
!$TCXOMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block) | ||
|
@@ -307,6 +375,100 @@ subroutine evp (dt) | |
enddo ! iblk | ||
!$TCXOMP END PARALLEL DO | ||
|
||
if (grid_system == 'CD') then | ||
|
||
!$TCXOMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block) | ||
do iblk = 1, nblocks | ||
|
||
!----------------------------------------------------------------- | ||
! more preparation for dynamics on N grid | ||
!----------------------------------------------------------------- | ||
|
||
this_block = get_block(blocks_ice(iblk),iblk) | ||
ilo = this_block%ilo | ||
ihi = this_block%ihi | ||
jlo = this_block%jlo | ||
jhi = this_block%jhi | ||
|
||
call dyn_prep2 (nx_block, ny_block, & | ||
ilo, ihi, jlo, jhi, & | ||
icellt(iblk), icelln(iblk), & | ||
indxti (:,iblk), indxtj (:,iblk), & | ||
indxni (:,iblk), indxnj (:,iblk), & | ||
aiN (:,:,iblk), nmass (:,:,iblk), & | ||
nmassdti (:,:,iblk), fcorN_blk (:,:,iblk), & | ||
nmask (:,:,iblk), & | ||
uocnN (:,:,iblk), vocnN (:,:,iblk), & | ||
strairxN (:,:,iblk), strairyN (:,:,iblk), & | ||
ss_tltxN (:,:,iblk), ss_tltyN (:,:,iblk), & | ||
icetmask (:,:,iblk), icenmask (:,:,iblk), & | ||
fmN (:,:,iblk), dt, & | ||
strtltxN (:,:,iblk), strtltyN (:,:,iblk), & | ||
strocnxN (:,:,iblk), strocnyN (:,:,iblk), & | ||
strintxN (:,:,iblk), strintyN (:,:,iblk), & | ||
taubxN (:,:,iblk), taubyN (:,:,iblk), & | ||
waterxN (:,:,iblk), wateryN (:,:,iblk), & | ||
forcexN (:,:,iblk), forceyN (:,:,iblk), & | ||
stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & | ||
stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & | ||
stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & | ||
stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & | ||
stress12_1(:,:,iblk), stress12_2(:,:,iblk), & | ||
stress12_3(:,:,iblk), stress12_4(:,:,iblk), & | ||
uvelN_init (:,:,iblk), vvelN_init (:,:,iblk), & | ||
uvelN (:,:,iblk), vvelN (:,:,iblk), & | ||
TbN (:,:,iblk)) | ||
|
||
enddo ! iblk | ||
!$TCXOMP END PARALLEL DO | ||
|
||
!$TCXOMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block) | ||
do iblk = 1, nblocks | ||
Comment on lines
+422
to
+426
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @dabail10 these two loops could be fusioned, right ? any reason why you added separate loops ? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This was fixed in PR #27 |
||
|
||
!----------------------------------------------------------------- | ||
! more preparation for dynamics on N grid | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. that should be: ! more preparation for dynamics on E grid |
||
!----------------------------------------------------------------- | ||
|
||
this_block = get_block(blocks_ice(iblk),iblk) | ||
ilo = this_block%ilo | ||
ihi = this_block%ihi | ||
jlo = this_block%jlo | ||
jhi = this_block%jhi | ||
|
||
call dyn_prep2 (nx_block, ny_block, & | ||
ilo, ihi, jlo, jhi, & | ||
icellt(iblk), icelle(iblk), & | ||
indxti (:,iblk), indxtj (:,iblk), & | ||
indxei (:,iblk), indxej (:,iblk), & | ||
aiE (:,:,iblk), emass (:,:,iblk), & | ||
emassdti (:,:,iblk), fcorE_blk (:,:,iblk), & | ||
emask (:,:,iblk), & | ||
uocnE (:,:,iblk), vocnE (:,:,iblk), & | ||
strairxE (:,:,iblk), strairyE (:,:,iblk), & | ||
ss_tltxE (:,:,iblk), ss_tltyE (:,:,iblk), & | ||
icetmask (:,:,iblk), icenmask (:,:,iblk), & | ||
fmE (:,:,iblk), dt, & | ||
strtltxE (:,:,iblk), strtltyE (:,:,iblk), & | ||
strocnxE (:,:,iblk), strocnyE (:,:,iblk), & | ||
strintxE (:,:,iblk), strintyE (:,:,iblk), & | ||
taubxE (:,:,iblk), taubyE (:,:,iblk), & | ||
waterxE (:,:,iblk), wateryE (:,:,iblk), & | ||
forcexE (:,:,iblk), forceyE (:,:,iblk), & | ||
stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & | ||
stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & | ||
stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & | ||
stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & | ||
stress12_1(:,:,iblk), stress12_2(:,:,iblk), & | ||
stress12_3(:,:,iblk), stress12_4(:,:,iblk), & | ||
uvelE_init (:,:,iblk), vvelE_init (:,:,iblk), & | ||
uvelE (:,:,iblk), vvelE (:,:,iblk), & | ||
TbE (:,:,iblk)) | ||
|
||
enddo ! iblk | ||
!$TCXOMP END PARALLEL DO | ||
|
||
endif ! grid_system | ||
|
||
call icepack_warnings_flush(nu_diag) | ||
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & | ||
file=__FILE__, line=__LINE__) | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -84,10 +84,22 @@ module ice_dyn_shared | |
real (kind=dbl_kind), allocatable, public :: & | ||
fcor_blk(:,:,:) ! Coriolis parameter (1/s) | ||
|
||
real (kind=dbl_kind), allocatable, public :: & | ||
fcorE_blk(:,:,:), & ! Coriolis parameter at E points (1/s) | ||
fcorN_blk(:,:,:) ! Coriolis parameter at N points (1/s) | ||
|
||
real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & | ||
uvel_init, & ! x-component of velocity (m/s), beginning of timestep | ||
vvel_init ! y-component of velocity (m/s), beginning of timestep | ||
|
||
real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & | ||
uvelN_init, & ! x-component of velocity (m/s), beginning of timestep | ||
vvelN_init ! y-component of velocity (m/s), beginning of timestep | ||
|
||
real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & | ||
uvelE_init, & ! x-component of velocity (m/s), beginning of timestep | ||
vvelE_init ! y-component of velocity (m/s), beginning of timestep | ||
|
||
! ice isotropic tensile strength parameter | ||
real (kind=dbl_kind), public :: & | ||
Ktens ! T=Ktens*P (tensile strength: see Konig and Holland, 2010) | ||
|
@@ -121,6 +133,16 @@ subroutine alloc_dyn_shared | |
stat=ierr) | ||
if (ierr/=0) call abort_ice('(alloc_dyn_shared): Out of memory') | ||
|
||
if (grid_system == 'CD') then | ||
allocate( & | ||
uvelE_init (nx_block,ny_block,max_blocks), & ! x-component of velocity (m/s), beginning of timestep | ||
vvelE_init (nx_block,ny_block,max_blocks), & ! y-component of velocity (m/s), beginning of timestep | ||
uvelN_init (nx_block,ny_block,max_blocks), & ! x-component of velocity (m/s), beginning of timestep | ||
vvelN_init (nx_block,ny_block,max_blocks), & ! y-component of velocity (m/s), beginning of timestep | ||
stat=ierr) | ||
if (ierr/=0) call abort_ice('(alloc_dyn_shared): Out of memory') | ||
endif | ||
|
||
end subroutine alloc_dyn_shared | ||
|
||
!======================================================================= | ||
|
@@ -138,7 +160,7 @@ subroutine init_dyn (dt) | |
stressm_1, stressm_2, stressm_3, stressm_4, & | ||
stress12_1, stress12_2, stress12_3, stress12_4 | ||
use ice_state, only: uvel, vvel, uvelE, vvelE, uvelN, vvelN, divu, shear | ||
use ice_grid, only: ULAT | ||
use ice_grid, only: ULAT, NLAT, ELAT | ||
|
||
real (kind=dbl_kind), intent(in) :: & | ||
dt ! time step | ||
|
@@ -161,6 +183,11 @@ subroutine init_dyn (dt) | |
|
||
allocate(fcor_blk(nx_block,ny_block,max_blocks)) | ||
|
||
if (grid_system == 'CD') then | ||
allocate(fcorE_blk(nx_block,ny_block,max_blocks)) | ||
allocate(fcorN_blk(nx_block,ny_block,max_blocks)) | ||
endif | ||
|
||
!$OMP PARALLEL DO PRIVATE(iblk,i,j) | ||
do iblk = 1, nblocks | ||
do j = 1, ny_block | ||
|
@@ -170,10 +197,10 @@ subroutine init_dyn (dt) | |
uvel(i,j,iblk) = c0 ! m/s | ||
vvel(i,j,iblk) = c0 ! m/s | ||
if (grid_system == 'CD') then ! extra velocity variables | ||
uvelE = c0 | ||
vvelE = c0 | ||
uvelN = c0 | ||
vvelN = c0 | ||
uvelE(i,j,iblk) = c0 | ||
vvelE(i,j,iblk) = c0 | ||
uvelN(i,j,iblk) = c0 | ||
vvelN(i,j,iblk) = c0 | ||
endif | ||
|
||
! strain rates | ||
|
@@ -191,6 +218,21 @@ subroutine init_dyn (dt) | |
fcor_blk(i,j,iblk) = c2*omega*sin(ULAT(i,j,iblk)) ! 1/s | ||
endif | ||
|
||
if (grid_system == 'CD') then | ||
|
||
if (trim(coriolis) == 'constant') then | ||
fcorE_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s | ||
fcorN_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s | ||
else if (trim(coriolis) == 'zero') then | ||
fcorE_blk(i,j,iblk) = 0.0 | ||
fcorN_blk(i,j,iblk) = 0.0 | ||
Comment on lines
+227
to
+228
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. these should be There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Good catch. This might change answers for zero coriolis tests. |
||
else | ||
fcorE_blk(i,j,iblk) = c2*omega*sin(ELAT(i,j,iblk)) ! 1/s | ||
fcorN_blk(i,j,iblk) = c2*omega*sin(NLAT(i,j,iblk)) ! 1/s | ||
endif | ||
|
||
endif | ||
|
||
! stress tensor, kg/s^2 | ||
stressp_1 (i,j,iblk) = c0 | ||
stressp_2 (i,j,iblk) = c0 | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I would prefer that we indent everything consistently. It makes code folding much more useful in editors, and makes reading the code easier (in my opinion)..