Skip to content
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.

Sign up for GitHub

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

(+) porous topography implementation #3

Merged
merged 9 commits into from
Nov 30, 2021
Merged

(+) porous topography implementation #3

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
8b5e10a
(+) porous topography implementation (#3)
sditkovsky Nov 17, 2021
bc82609
small fixes to porous topo code
sditkovsky Nov 19, 2021
e717e46
preserve dimensions in porous topo interpolation
sditkovsky Nov 24, 2021
04d932e
Remove Z_to_m rescaling in porous module
sditkovsky Nov 24, 2021
a98c3ef
small bug fix
sditkovsky Nov 24, 2021
66e1493
small bug fix 2
sditkovsky Nov 24, 2021
62786c4
fix merge issue
sditkovsky Nov 24, 2021
ec84be5
change conversion H_to_MKS to H_to_m for u/veffA
sditkovsky Nov 29, 2021
3f38239
Merge branch 'dev/gfdl' into porous_topo_sjd
sditkovsky Nov 29, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/core/MOM_dynamics_split_RK2.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -398,7 +398,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
do j=G%jsdB,G%jedB ; do i=G%isd,G%ied ; vp(i,j,k) = 0.0 ; enddo ; enddo
do j=G%jsd,G%jed ; do i=G%isd,G%ied ; hp(i,j,k) = h(i,j,k) ; enddo ; enddo
enddo
ueffA(:,:,:) = 0; veffA(:,:,:) = 0
if (CS%id_ueffA > 0) ueffA(:,:,:) = 0
if (CS%id_veffA > 0) veffA(:,:,:) = 0

! Update CFL truncation value as function of time
call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp)
Expand Down
3 changes: 2 additions & 1 deletion src/core/MOM_dynamics_unsplit.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -244,7 +244,8 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, &
h_av(:,:,:) = 0; hp(:,:,:) = 0
up(:,:,:) = 0; upp(:,:,:) = 0
vp(:,:,:) = 0; vpp(:,:,:) = 0
ueffA(:,:,:) = 0; veffA(:,:,:) = 0
if (CS%id_ueffA > 0) ueffA(:,:,:) = 0
if (CS%id_veffA > 0) veffA(:,:,:) = 0

dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
if (dyn_p_surf) then
Expand Down
3 changes: 2 additions & 1 deletion src/core/MOM_dynamics_unsplit_RK2.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -254,7 +254,8 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
h_av(:,:,:) = 0; hp(:,:,:) = 0
up(:,:,:) = 0
vp(:,:,:) = 0
ueffA(:,:,:) = 0; veffA(:,:,:) = 0
if (CS%id_ueffA > 0) ueffA(:,:,:) = 0
if (CS%id_veffA > 0) veffA(:,:,:) = 0

dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
if (dyn_p_surf) then
Expand Down
12 changes: 6 additions & 6 deletions src/core/MOM_grid.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,14 +113,14 @@ module MOM_grid
areaCv !< The areas of the v-grid cells [L2 ~> m2].

real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: &
porous_DminU, & !< minimum topographic height of U-face [m]
porous_DmaxU, & !< maximum topographic height of U-face [m]
porous_DavgU !< average topographic height of U-face [m]
porous_DminU, & !< minimum topographic height of U-face [Z ~> m]
porous_DmaxU, & !< maximum topographic height of U-face [Z ~> m]
porous_DavgU !< average topographic height of U-face [Z ~> m]

real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: &
porous_DminV, & !< minimum topographic height of V-face [m]
porous_DmaxV, & !< maximum topographic height of V-face [m]
porous_DavgV !< average topographic height of V-face [m]
porous_DminV, & !< minimum topographic height of V-face [Z ~> m]
porous_DmaxV, & !< maximum topographic height of V-face [Z ~> m]
porous_DavgV !< average topographic height of V-face [Z ~> m]

real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
mask2dBu, & !< 0 for boundary points and 1 for ocean points on the q grid [nondim].
Expand Down
138 changes: 72 additions & 66 deletions src/core/MOM_porous_barriers.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
!> Function for calculating curve fit for porous topography.
!> Module for calculating curve fit for porous topography.
!written by sjd
module MOM_porous_barriers

Expand All @@ -24,6 +24,7 @@ module MOM_porous_barriers

contains

!> subroutine to assign cell face areas and layer widths for porous topography
subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m)
!eta_bt, halo_size, eta_to_m not currently used
!variables needed to call find_eta
Expand Down Expand Up @@ -57,69 +58,74 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m)
IsdB = G%IsdB; IedB = G%IedB; JsdB = G%JsdB; JedB = G%JedB

!eta is zero at surface and decreases downward
!all calculations are done in [m]

nk = SZK_(G)

!currently no treatment for using optional find_eta arguments if present
call find_eta(h, tv, G, GV, US, eta)

do I=IsdB,IedB; do j=jsd,jed
if (G%porous_DavgU(I,j) < 0.) then
do K = nk+1,1,-1
eta_s = max(US%Z_to_m*eta(I,j,K), US%Z_to_m*eta(I+1,j,K)) !take shallower layer height
!eta_s = 0.5 * (US%Z_to_m*eta(I,j,K) + US%Z_to_m*eta(I+1,j,K)) !take arithmetic mean
if (eta_s <= G%porous_DminU(I,j)) then
pbv%por_layer_widthU(I,j,K) = 0.0
A_layer_prev = 0.0
if (K < nk+1) then
pbv%por_face_areaU(I,j,k) = 0.0; endif
else
call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j),&
eta_s, w_layer, A_layer)
pbv%por_layer_widthU(I,j,K) = w_layer
if (k <= nk) then
if ((eta_s - eta_prev) > 0.0) then
pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev)/&
do j=jsd,jed; do I=IsdB,IedB
if (G%porous_DavgU(I,j) < 0.) then
do K = nk+1,1,-1
eta_s = max(US%Z_to_m*eta(I,j,K), US%Z_to_m*eta(I+1,j,K)) !take shallower layer height
!eta_s = 0.5 * (US%Z_to_m*eta(I,j,K) + US%Z_to_m*eta(I+1,j,K)) !take arithmetic mean
if (eta_s <= G%porous_DminU(I,j)) then
pbv%por_layer_widthU(I,j,K) = 0.0
A_layer_prev = 0.0
if (K < nk+1) then
pbv%por_face_areaU(I,j,k) = 0.0; endif
else
call calc_por_layer(US%Z_to_m*(G%porous_DminU(I,j)-G%Z_ref), &
US%Z_to_m*(G%porous_DmaxU(I,j)-G%Z_ref), &
US%Z_to_m*(G%porous_DavgU(I,j)-G%Z_ref), eta_s, w_layer, A_layer)
pbv%por_layer_widthU(I,j,K) = w_layer
if (k <= nk) then
if ((eta_s - eta_prev) > 0.0) then
pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev)/&
(eta_s-eta_prev)
else
pbv%por_face_areaU(I,j,k) = 0.0; endif
endif
eta_prev = eta_s
A_layer_prev = A_layer
endif; enddo
endif; enddo; enddo
else
pbv%por_face_areaU(I,j,k) = 0.0; endif
endif
eta_prev = eta_s
A_layer_prev = A_layer
endif
enddo
endif
enddo; enddo

do i=isd,ied; do J=JsdB,JedB
if (G%porous_DavgV(i,J) < 0.) then
do K = nk+1,1,-1
eta_s = max(US%Z_to_m*eta(i,J,K), US%Z_to_m*eta(i,J+1,K)) !take shallower layer height
!eta_s = 0.5 * (US%Z_to_m*eta(i,J,K) + US%Z_to_m*eta(i,J+1,K)) !take arithmetic mean
if (eta_s <= G%porous_DminV(i,J)) then
pbv%por_layer_widthV(i,J,K) = 0.0
A_layer_prev = 0.0
if (K < nk+1) then
pbv%por_face_areaV(i,J,k) = 0.0; endif
else
call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J),&
eta_s, w_layer, A_layer)
pbv%por_layer_widthV(i,J,K) = w_layer
if (k <= nk) then
if ((eta_s - eta_prev) > 0.0) then
pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev)/&
if (G%porous_DavgV(i,J) < 0.) then
do K = nk+1,1,-1
eta_s = max(US%Z_to_m*eta(i,J,K), US%Z_to_m*eta(i,J+1,K)) !take shallower layer height
!eta_s = 0.5 * (US%Z_to_m*eta(i,J,K) + US%Z_to_m*eta(i,J+1,K)) !take arithmetic mean
if (eta_s <= G%porous_DminV(i,J)) then
pbv%por_layer_widthV(i,J,K) = 0.0
A_layer_prev = 0.0
if (K < nk+1) then
pbv%por_face_areaV(i,J,k) = 0.0; endif
else
call calc_por_layer(US%Z_to_m*(G%porous_DminV(i,J)-G%Z_ref), &
US%Z_to_m*(G%porous_DmaxV(i,J)-G%Z_ref), &
US%Z_to_m*(G%porous_DavgV(i,J)-G%Z_ref), eta_s, w_layer, A_layer)
pbv%por_layer_widthV(i,J,K) = w_layer
if (k <= nk) then
if ((eta_s - eta_prev) > 0.0) then
pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev)/&
(eta_s-eta_prev)
else
pbv%por_face_areaU(I,j,k) = 0.0; endif
endif
eta_prev = eta_s
A_layer_prev = A_layer
endif; enddo
endif; enddo; enddo
else
pbv%por_face_areaU(I,j,k) = 0.0; endif
endif
eta_prev = eta_s
A_layer_prev = A_layer
endif
enddo
endif
enddo; enddo

end subroutine por_widths

!> subroutine to calculate the profile fit for a single layer in a column
subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer)
!subroutine to calculate the profile fit for a layer

real, intent(in) :: D_min !< minimum topographic height [m]
real, intent(in) :: D_max !< maximum topographic height [m]
Expand All @@ -140,24 +146,24 @@ subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer)
zeta = (eta_layer - D_min)/(D_max - D_min)

if (eta_layer <= D_min) then
w_layer = 0.0
A_layer = 0.0
w_layer = 0.0
A_layer = 0.0
elseif (eta_layer >= D_max) then
w_layer = 1.0
A_layer = eta_layer - D_avg
w_layer = 1.0
A_layer = eta_layer - D_avg
else
if (m < 0.5) then
psi = zeta**(1./a)
psi_int = (1.-m)*zeta**(1./(1.-m))
elseif (m == 0.5) then
psi = zeta
psi_int = 0.5*zeta*zeta
else
psi = 1. - (1. - zeta)**a
psi_int = zeta - m + m*((1-zeta)**(1/m))
endif
w_layer = psi
A_layer = (D_max - D_min)*psi_int
if (m < 0.5) then
psi = zeta**(1./a)
psi_int = (1.-m)*zeta**(1./(1.-m))
elseif (m == 0.5) then
psi = zeta
psi_int = 0.5*zeta*zeta
else
psi = 1. - (1. - zeta)**a
psi_int = zeta - m + m*((1-zeta)**(1/m))
endif
w_layer = psi
A_layer = (D_max - D_min)*psi_int
endif


Expand Down
24 changes: 12 additions & 12 deletions src/core/MOM_transcribe_grid.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,9 +71,9 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
oG%dyCu(I,j) = dG%dyCu(I+ido,j+jdo)
oG%dy_Cu(I,j) = dG%dy_Cu(I+ido,j+jdo)

oG%porous_DminU(I,j) = dG%porous_DminU(I+ido,j+jdo)
oG%porous_DmaxU(I,j) = dG%porous_DmaxU(I+ido,j+jdo)
oG%porous_DavgU(I,j) = dG%porous_DavgU(I+ido,j+jdo)
oG%porous_DminU(I,j) = dG%porous_DminU(I+ido,j+jdo) - oG%Z_ref
oG%porous_DmaxU(I,j) = dG%porous_DmaxU(I+ido,j+jdo) - oG%Z_ref
oG%porous_DavgU(I,j) = dG%porous_DavgU(I+ido,j+jdo) - oG%Z_ref

oG%mask2dCu(I,j) = dG%mask2dCu(I+ido,j+jdo)
oG%areaCu(I,j) = dG%areaCu(I+ido,j+jdo)
Expand All @@ -87,9 +87,9 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
oG%dyCv(i,J) = dG%dyCv(i+ido,J+jdo)
oG%dx_Cv(i,J) = dG%dx_Cv(i+ido,J+jdo)

oG%porous_DminV(i,J) = dG%porous_DminV(i+ido,J+jdo)
oG%porous_DmaxV(i,J) = dG%porous_DmaxV(i+ido,J+jdo)
oG%porous_DavgV(i,J) = dG%porous_DavgV(i+ido,J+jdo)
oG%porous_DminV(i,J) = dG%porous_DminV(i+ido,J+jdo) - oG%Z_ref
oG%porous_DmaxV(i,J) = dG%porous_DmaxV(i+ido,J+jdo) - oG%Z_ref
oG%porous_DavgV(i,J) = dG%porous_DavgV(i+ido,J+jdo) - oG%Z_ref

oG%mask2dCv(i,J) = dG%mask2dCv(i+ido,J+jdo)
oG%areaCv(i,J) = dG%areaCv(i+ido,J+jdo)
Expand Down Expand Up @@ -224,9 +224,9 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
dG%dyCu(I,j) = oG%dyCu(I+ido,j+jdo)
dG%dy_Cu(I,j) = oG%dy_Cu(I+ido,j+jdo)

dG%porous_DminU(I,j) = oG%porous_DminU(I+ido,j+jdo)
dG%porous_DmaxU(I,j) = oG%porous_DmaxU(I+ido,j+jdo)
dG%porous_DavgU(I,j) = oG%porous_DavgU(I+ido,j+jdo)
dG%porous_DminU(I,j) = oG%porous_DminU(I+ido,j+jdo) + oG%Z_ref
dG%porous_DmaxU(I,j) = oG%porous_DmaxU(I+ido,j+jdo) + oG%Z_ref
dG%porous_DavgU(I,j) = oG%porous_DavgU(I+ido,j+jdo) + oG%Z_ref

dG%mask2dCu(I,j) = oG%mask2dCu(I+ido,j+jdo)
dG%areaCu(I,j) = oG%areaCu(I+ido,j+jdo)
Expand All @@ -240,9 +240,9 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
dG%dyCv(i,J) = oG%dyCv(i+ido,J+jdo)
dG%dx_Cv(i,J) = oG%dx_Cv(i+ido,J+jdo)

dG%porous_DminV(i,J) = oG%porous_DminU(i+ido,J+jdo)
dG%porous_DmaxV(i,J) = oG%porous_DmaxU(i+ido,J+jdo)
dG%porous_DavgV(i,J) = oG%porous_DavgU(i+ido,J+jdo)
dG%porous_DminV(i,J) = oG%porous_DminU(i+ido,J+jdo) + oG%Z_ref
dG%porous_DmaxV(i,J) = oG%porous_DmaxU(i+ido,J+jdo) + oG%Z_ref
dG%porous_DavgV(i,J) = oG%porous_DavgU(i+ido,J+jdo) + oG%Z_ref

dG%mask2dCv(i,J) = oG%mask2dCv(i+ido,J+jdo)
dG%areaCv(i,J) = oG%areaCv(i+ido,J+jdo)
Expand Down
12 changes: 6 additions & 6 deletions src/framework/MOM_dyn_horgrid.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,14 +110,14 @@ module MOM_dyn_horgrid
areaCv !< The areas of the v-grid cells [L2 ~> m2].

real, allocatable, dimension(:,:) :: &
porous_DminU, & !< minimum topographic height of U-face [m]
porous_DmaxU, & !< maximum topographic height of U-face [m]
porous_DavgU !< average topographic height of U-face [m]
porous_DminU, & !< minimum topographic height of U-face [Z ~> m]
porous_DmaxU, & !< maximum topographic height of U-face [Z ~> m]
porous_DavgU !< average topographic height of U-face [Z ~> m]

real, allocatable, dimension(:,:) :: &
porous_DminV, & !< minimum topographic height of V-face [m]
porous_DmaxV, & !< maximum topographic height of V-face [m]
porous_DavgV !< average topographic height of V-face [m]
porous_DminV, & !< minimum topographic height of V-face [Z ~> m]
porous_DmaxV, & !< maximum topographic height of V-face [Z ~> m]
porous_DavgV !< average topographic height of V-face [Z ~> m]

real, allocatable, dimension(:,:) :: &
mask2dBu, & !< 0 for boundary points and 1 for ocean points on the q grid [nondim].
Expand Down
22 changes: 12 additions & 10 deletions src/initialization/MOM_shared_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -856,6 +856,7 @@ subroutine reset_face_lengths_list(G, param_file, US)
Dmin_v, Dmax_v, Davg_v
real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim]
real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim]
real :: m_to_Z ! A unit conversion factor [L ~> m]
real :: lat, lon ! The latitude and longitude of a point.
real :: len_lon ! The periodic range of longitudes, usually 360 degrees.
real :: len_lat ! The range of latitudes, usually 180 degrees.
Expand All @@ -878,6 +879,7 @@ subroutine reset_face_lengths_list(G, param_file, US)
call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L
L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m
m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z

call get_param(param_file, mdl, "CHANNEL_LIST_FILE", chan_file, &
"The file from which the list of narrowed channels is read.", &
Expand Down Expand Up @@ -971,9 +973,9 @@ subroutine reset_face_lengths_list(G, param_file, US)
if (found_u) then
u_pt = u_pt + 1
if (found_u_por .eqv. .false.) then
read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt)
read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt)
elseif (found_u_por) then
read(line(isu_por+12:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt), &
read(line(isu_por+12:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt), &
Dmin_u(u_pt), Dmax_u(u_pt), Davg_u(u_pt)
endif
u_line_no(u_pt) = ln
Expand Down Expand Up @@ -1008,9 +1010,9 @@ subroutine reset_face_lengths_list(G, param_file, US)
elseif (found_v) then
v_pt = v_pt + 1
if (found_v_por .eqv. .false.) then
read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt)
read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt)
elseif (found_v_por) then
read(line(isv+12:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt), &
read(line(isv+12:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt), &
Dmin_v(v_pt), Dmax_v(v_pt), Davg_v(v_pt)
endif
v_line_no(v_pt) = ln
Expand Down Expand Up @@ -1059,9 +1061,9 @@ subroutine reset_face_lengths_list(G, param_file, US)
((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) ) then

G%dy_Cu(I,j) = G%mask2dCu(I,j) * m_to_L*min(L_to_m*G%dyCu(I,j), max(u_width(npt), 0.0))
G%porous_DminU(I,j) = Dmin_u(npt)
G%porous_DmaxU(I,j) = Dmax_u(npt)
G%porous_DavgU(I,j) = Davg_u(npt)
G%porous_DminU(I,j) = m_to_Z*Dmin_u(npt)
G%porous_DmaxU(I,j) = m_to_Z*Dmax_u(npt)
G%porous_DavgU(I,j) = m_to_Z*Davg_u(npt)

if (j>=G%jsc .and. j<=G%jec .and. I>=G%isc .and. I<=G%iec) then ! Limit messages/checking to compute domain
if ( G%mask2dCu(I,j) == 0.0 ) then
Expand Down Expand Up @@ -1096,9 +1098,9 @@ subroutine reset_face_lengths_list(G, param_file, US)
((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. &
((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) ) then
G%dx_Cv(i,J) = G%mask2dCv(i,J) * m_to_L*min(L_to_m*G%dxCv(i,J), max(v_width(npt), 0.0))
G%porous_DminV(i,J) = Dmin_v(npt)
G%porous_DmaxV(i,J) = Dmax_v(npt)
G%porous_DavgV(i,J) = Davg_v(npt)
G%porous_DminV(i,J) = m_to_Z*Dmin_v(npt)
G%porous_DmaxV(i,J) = m_to_Z*Dmax_v(npt)
G%porous_DavgV(i,J) = m_to_Z*Davg_v(npt)
if (i>=G%isc .and. i<=G%iec .and. J>=G%jsc .and. J<=G%jec) then ! Limit messages/checking to compute domain
if ( G%mask2dCv(i,J) == 0.0 ) then
write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",&
Expand Down