diff --git a/ext/AGRIF/AGRIF_FILES/modbcfunction.F90 b/ext/AGRIF/AGRIF_FILES/modbcfunction.F90
index e6add0329e709fbd64311d119d638024a7ad141d..ae7de45aa4f67ad6544725c894575005c8d1f88c 100644
--- a/ext/AGRIF/AGRIF_FILES/modbcfunction.F90
+++ b/ext/AGRIF/AGRIF_FILES/modbcfunction.F90
@@ -447,14 +447,14 @@ subroutine Agrif_Find_Nearest ( tabvarsindic, fineloc, parentloc )
type(Agrif_Variable), pointer :: parent_var
type(Agrif_Variable), pointer :: child_var
integer :: i
- integer, dimension(6) :: nb_child !< Number of cells on the child grid
- integer, dimension(6) :: ub_child !< Upper bound on the child grid
- integer, dimension(6) :: lb_child !< Lower bound on the child grid
- integer, dimension(6) :: lb_parent !< Lower bound on the parent grid
- real, dimension(6) :: s_child !< Child grid position (s_root = 0)
- real, dimension(6) :: s_parent !< Parent grid position (s_root = 0)
- real, dimension(6) :: ds_child !< Child grid dx (ds_root = 1)
- real, dimension(6) :: ds_parent !< Parent grid dx (ds_root = 1)
+ integer, dimension(6) :: nb_child !< Number of cells on the child grid
+ integer, dimension(6) :: ub_child !< Upper bound on the child grid
+ integer, dimension(6) :: lb_child !< Lower bound on the child grid
+ integer, dimension(6) :: lb_parent !< Lower bound on the parent grid
+ real(kind=8), dimension(6) :: s_child !< Child grid position (s_root = 0)
+ real(kind=8), dimension(6) :: s_parent !< Parent grid position (s_root = 0)
+ real(kind=8), dimension(6) :: ds_child !< Child grid dx (ds_root = 1)
+ real(kind=8), dimension(6) :: ds_parent !< Parent grid dx (ds_root = 1)
integer :: nbdim !< Number of dimensions
real, dimension(6) :: xfineloc
!
diff --git a/ext/AGRIF/AGRIF_FILES/modupdate.F90 b/ext/AGRIF/AGRIF_FILES/modupdate.F90
index 3bdb5b1b32c32a04011776637e240d6645e4d5c5..c647fcb7930460be33cc07e7cb3b51fdb5f3e63d 100644
--- a/ext/AGRIF/AGRIF_FILES/modupdate.F90
+++ b/ext/AGRIF/AGRIF_FILES/modupdate.F90
@@ -1904,8 +1904,8 @@ subroutine Agrif_Update_1D_Recursive ( type_update, &
integer, intent(in) :: type_update !< Type of update (copy or average)
integer, intent(in) :: indmin, indmax
integer, intent(in) :: lb_child, ub_child
- real(kind=8), intent(in) :: s_child, s_parent
- real(kind=8), intent(in) :: ds_child, ds_parent
+ real(kind=8), intent(in) :: s_child, s_parent
+ real(kind=8), intent(in) :: ds_child, ds_parent
real, dimension(indmin:indmax), intent(out) :: tempP
real, dimension(lb_child:ub_child), intent(in) :: tempC
!---------------------------------------------------------------------------------------------------
@@ -1933,11 +1933,11 @@ subroutine Agrif_Update_2D_Recursive ( type_update, &
s_child, s_parent, &
ds_child, ds_parent )
!---------------------------------------------------------------------------------------------------
- integer, dimension(2), intent(in) :: type_update !< Type of update (copy or average)
- integer, dimension(2), intent(in) :: indmin, indmax
- integer, dimension(2), intent(in) :: lb_child, ub_child
- real(kind=8), dimension(2), intent(in) :: s_child, s_parent
- real(kind=8), dimension(2), intent(in) :: ds_child, ds_parent
+ integer, dimension(2), intent(in) :: type_update !< Type of update (copy or average)
+ integer, dimension(2), intent(in) :: indmin, indmax
+ integer, dimension(2), intent(in) :: lb_child, ub_child
+ real(kind=8),dimension(2), intent(in) :: s_child, s_parent
+ real(kind=8),dimension(2), intent(in) :: ds_child, ds_parent
real, dimension( &
indmin(1):indmax(1), &
indmin(2):indmax(2)), intent(out) :: tempP
@@ -2047,11 +2047,11 @@ subroutine Agrif_Update_2D_Recursive_ok ( type_update, &
lb_child, ub_child, &
s_child, s_parent, ds_child, ds_parent )
!---------------------------------------------------------------------------------------------------
- INTEGER, DIMENSION(2), intent(in) :: type_update !< Type of update (copy or average)
- INTEGER, DIMENSION(2), intent(in) :: indmin, indmax
- INTEGER, DIMENSION(2), intent(in) :: lb_child, ub_child
- REAL, DIMENSION(2), intent(in) :: s_child, s_parent
- REAL, DIMENSION(2), intent(in) :: ds_child, ds_parent
+ INTEGER, DIMENSION(2), intent(in) :: type_update !< Type of update (copy or average)
+ INTEGER, DIMENSION(2), intent(in) :: indmin, indmax
+ INTEGER, DIMENSION(2), intent(in) :: lb_child, ub_child
+ REAL(kind=8),DIMENSION(2), intent(in) :: s_child, s_parent
+ REAL(kind=8),DIMENSION(2), intent(in) :: ds_child, ds_parent
REAL, DIMENSION( &
indmin(1):indmax(1), &
indmin(2):indmax(2)), intent(out) :: tempP
@@ -2101,11 +2101,11 @@ subroutine Agrif_Update_3D_Recursive ( type_update, &
s_child, s_parent, &
ds_child, ds_parent )
!---------------------------------------------------------------------------------------------------
- integer, dimension(3), intent(in) :: type_update !< Type of update (copy or average)
- integer, dimension(3), intent(in) :: indmin, indmax
- integer, dimension(3), intent(in) :: lb_child, ub_child
- real, dimension(3), intent(in) :: s_child, s_parent
- real, dimension(3), intent(in) :: ds_child, ds_parent
+ integer, dimension(3), intent(in) :: type_update !< Type of update (copy or average)
+ integer, dimension(3), intent(in) :: indmin, indmax
+ integer, dimension(3), intent(in) :: lb_child, ub_child
+ real(kind=8),dimension(3), intent(in) :: s_child, s_parent
+ real(kind=8),dimension(3), intent(in) :: ds_child, ds_parent
real, dimension( &
indmin(1):indmax(1), &
indmin(2):indmax(2), &
@@ -2213,11 +2213,11 @@ subroutine Agrif_Update_4D_Recursive ( type_update, &
s_child, s_parent, &
ds_child, ds_parent )
!---------------------------------------------------------------------------------------------------
- integer, dimension(4), intent(in) :: type_update !< Type of update (copy or average)
- integer, dimension(4), intent(in) :: indmin, indmax
- integer, dimension(4), intent(in) :: lb_child, ub_child
- real, dimension(4), intent(in) :: s_child, s_parent
- real, dimension(4), intent(in) :: ds_child, ds_parent
+ integer, dimension(4), intent(in) :: type_update !< Type of update (copy or average)
+ integer, dimension(4), intent(in) :: indmin, indmax
+ integer, dimension(4), intent(in) :: lb_child, ub_child
+ real(kind=8),dimension(4), intent(in) :: s_child, s_parent
+ real(kind=8),dimension(4), intent(in) :: ds_child, ds_parent
real, dimension( &
indmin(1):indmax(1), &
indmin(2):indmax(2), &
@@ -2285,11 +2285,11 @@ subroutine Agrif_Update_5D_Recursive ( type_update, &
s_child, s_parent, &
ds_child, ds_parent )
!---------------------------------------------------------------------------------------------------
- integer, dimension(5), intent(in) :: type_update !< Type of update (copy or average)
- integer, dimension(5), intent(in) :: indmin, indmax
- integer, dimension(5), intent(in) :: lb_child, ub_child
- real, dimension(5), intent(in) :: s_child, s_parent
- real, dimension(5), intent(in) :: ds_child, ds_parent
+ integer, dimension(5), intent(in) :: type_update !< Type of update (copy or average)
+ integer, dimension(5), intent(in) :: indmin, indmax
+ integer, dimension(5), intent(in) :: lb_child, ub_child
+ real(kind=8),dimension(5), intent(in) :: s_child, s_parent
+ real(kind=8),dimension(5), intent(in) :: ds_child, ds_parent
real, dimension( &
indmin(1):indmax(1), &
indmin(2):indmax(2), &
@@ -2364,11 +2364,11 @@ subroutine Agrif_Update_6D_Recursive ( type_update, &
s_child, s_parent, &
ds_child, ds_parent )
!---------------------------------------------------------------------------------------------------
- integer, dimension(6), intent(in) :: type_update !< Type of update (copy or average)
- integer, dimension(6), intent(in) :: indmin, indmax
- integer, dimension(6), intent(in) :: lb_child, ub_child
- real, dimension(6), intent(in) :: s_child, s_parent
- real, dimension(6), intent(in) :: ds_child, ds_parent
+ integer, dimension(6), intent(in) :: type_update !< Type of update (copy or average)
+ integer, dimension(6), intent(in) :: indmin, indmax
+ integer, dimension(6), intent(in) :: lb_child, ub_child
+ real(kind=8),dimension(6), intent(in) :: s_child, s_parent
+ real(kind=8),dimension(6), intent(in) :: ds_child, ds_parent
real, dimension( &
indmin(1):indmax(1), &
indmin(2):indmax(2), &
@@ -2454,8 +2454,8 @@ subroutine Agrif_UpdateBase ( type_update, &
integer, intent(in) :: lb_child, ub_child
real, dimension(indmin:indmax), intent(out):: parent_tab
real, dimension(lb_child:ub_child), intent(in) :: child_tab
- real, intent(in) :: s_parent, s_child
- real, intent(in) :: ds_parent, ds_child
+ real(kind=8), intent(in) :: s_parent, s_child
+ real(kind=8), intent(in) :: ds_parent, ds_child
!---------------------------------------------------------------------------------------------------
integer :: np ! Length of parent array
integer :: nc ! Length of child array
diff --git a/ext/AGRIF/AGRIF_FILES/modupdatebasic.F90 b/ext/AGRIF/AGRIF_FILES/modupdatebasic.F90
index 900f0386f97c1dc62b574a4172085baf9759eb9d..76367bcfe79fbba0a5831ad285fd57b26f351332 100644
--- a/ext/AGRIF/AGRIF_FILES/modupdatebasic.F90
+++ b/ext/AGRIF/AGRIF_FILES/modupdatebasic.F90
@@ -231,16 +231,15 @@ subroutine Agrif_basicupdate_max1d ( x, y, np, nc, s_parent, s_child, ds_parent,
REAL, DIMENSION(np), intent(out) :: x
REAL, DIMENSION(nc), intent(in) :: y
INTEGER, intent(in) :: np,nc
- REAL, intent(in) :: s_parent, s_child
- REAL, intent(in) :: ds_parent, ds_child
+ REAL(kind=8), intent(in) :: s_parent, s_child
+ REAL(kind=8), intent(in) :: ds_parent, ds_child
!
INTEGER :: i, ii, locind_child_left, coeffraf
- REAL :: xpos, invcoeffraf
+ REAL(kind=8) :: xpos
INTEGER :: nbnonnuls
INTEGER :: diffmod
!
coeffraf = nint(ds_parent/ds_child)
- invcoeffraf = 1./coeffraf
!
if (coeffraf == 1) then
locind_child_left = 1 + nint((s_parent - s_child)/ds_child)
diff --git a/ext/PPR/src/bfun1d.h90 b/ext/PPR/src/bfun1d.h90
index 279634678bb4805d2e674e5479a3e09bbaa0b257..1a40a981c6d3a4d69449580957363d1507761767 100644
--- a/ext/PPR/src/bfun1d.h90
+++ b/ext/PPR/src/bfun1d.h90
@@ -59,46 +59,46 @@
!------------------------------------ -1th-order basis !
select case (ndof)
case (+1)
- bfun(1) = sval**1 / 1.e0
+ bfun(1) = sval**1 / 1.e0_8
case (+2)
- bfun(1) = sval**1 / 1.e0
- bfun(2) = sval**2 / 2.e0
+ bfun(1) = sval**1 / 1.e0_8
+ bfun(2) = sval**2 / 2.e0_8
case (+3)
- bfun(1) = sval**1 / 1.e0
- bfun(2) = sval**2 / 2.e0
- bfun(3) = sval**3 / 3.e0
+ bfun(1) = sval**1 / 1.e0_8
+ bfun(2) = sval**2 / 2.e0_8
+ bfun(3) = sval**3 / 3.e0_8
case (+4)
- bfun(1) = sval**1 / 1.e0
- bfun(2) = sval**2 / 2.e0
- bfun(3) = sval**3 / 3.e0
- bfun(4) = sval**4 / 4.e0
+ bfun(1) = sval**1 / 1.e0_8
+ bfun(2) = sval**2 / 2.e0_8
+ bfun(3) = sval**3 / 3.e0_8
+ bfun(4) = sval**4 / 4.e0_8
case (+5)
- bfun(1) = sval**1 / 1.e0
- bfun(2) = sval**2 / 2.e0
- bfun(3) = sval**3 / 3.e0
- bfun(4) = sval**4 / 4.e0
- bfun(5) = sval**5 / 5.e0
+ bfun(1) = sval**1 / 1.e0_8
+ bfun(2) = sval**2 / 2.e0_8
+ bfun(3) = sval**3 / 3.e0_8
+ bfun(4) = sval**4 / 4.e0_8
+ bfun(5) = sval**5 / 5.e0_8
case (+6)
- bfun(1) = sval**1 / 1.e0
- bfun(2) = sval**2 / 2.e0
- bfun(3) = sval**3 / 3.e0
- bfun(4) = sval**4 / 4.e0
- bfun(5) = sval**5 / 5.e0
- bfun(6) = sval**6 / 6.e0
+ bfun(1) = sval**1 / 1.e0_8
+ bfun(2) = sval**2 / 2.e0_8
+ bfun(3) = sval**3 / 3.e0_8
+ bfun(4) = sval**4 / 4.e0_8
+ bfun(5) = sval**5 / 5.e0_8
+ bfun(6) = sval**6 / 6.e0_8
case (+7)
- bfun(1) = sval**1 / 1.e0
- bfun(2) = sval**2 / 2.e0
- bfun(3) = sval**3 / 3.e0
- bfun(4) = sval**4 / 4.e0
- bfun(5) = sval**5 / 5.e0
- bfun(6) = sval**6 / 6.e0
- bfun(7) = sval**7 / 7.e0
+ bfun(1) = sval**1 / 1.e0_8
+ bfun(2) = sval**2 / 2.e0_8
+ bfun(3) = sval**3 / 3.e0_8
+ bfun(4) = sval**4 / 4.e0_8
+ bfun(5) = sval**5 / 5.e0_8
+ bfun(6) = sval**6 / 6.e0_8
+ bfun(7) = sval**7 / 7.e0_8
end select
@@ -106,46 +106,46 @@
!------------------------------------ +0th-order basis !
select case (ndof)
case (+1)
- bfun(1) = 1.e0
+ bfun(1) = 1.e0_8
case (+2)
- bfun(1) = 1.e0
- bfun(2) = sval**1 * 1.e0
+ bfun(1) = 1.e0_8
+ bfun(2) = sval**1 * 1.e0_8
case (+3)
- bfun(1) = 1.e0
- bfun(2) = sval**1 * 1.e0
- bfun(3) = sval**2 * 1.e0
+ bfun(1) = 1.e0_8
+ bfun(2) = sval**1 * 1.e0_8
+ bfun(3) = sval**2 * 1.e0_8
case (+4)
- bfun(1) = 1.e0
- bfun(2) = sval**1 * 1.e0
- bfun(3) = sval**2 * 1.e0
- bfun(4) = sval**3 * 1.e0
+ bfun(1) = 1.e0_8
+ bfun(2) = sval**1 * 1.e0_8
+ bfun(3) = sval**2 * 1.e0_8
+ bfun(4) = sval**3 * 1.e0_8
case (+5)
- bfun(1) = 1.e0
- bfun(2) = sval**1 * 1.e0
- bfun(3) = sval**2 * 1.e0
- bfun(4) = sval**3 * 1.e0
- bfun(5) = sval**4 * 1.e0
+ bfun(1) = 1.e0_8
+ bfun(2) = sval**1 * 1.e0_8
+ bfun(3) = sval**2 * 1.e0_8
+ bfun(4) = sval**3 * 1.e0_8
+ bfun(5) = sval**4 * 1.e0_8
case (+6)
- bfun(1) = 1.e0
- bfun(2) = sval**1 * 1.e0
- bfun(3) = sval**2 * 1.e0
- bfun(4) = sval**3 * 1.e0
- bfun(5) = sval**4 * 1.e0
- bfun(6) = sval**5 * 1.e0
+ bfun(1) = 1.e0_8
+ bfun(2) = sval**1 * 1.e0_8
+ bfun(3) = sval**2 * 1.e0_8
+ bfun(4) = sval**3 * 1.e0_8
+ bfun(5) = sval**4 * 1.e0_8
+ bfun(6) = sval**5 * 1.e0_8
case (+7)
- bfun(1) = 1.e0
- bfun(2) = sval**1 * 1.e0
- bfun(3) = sval**2 * 1.e0
- bfun(4) = sval**3 * 1.e0
- bfun(5) = sval**4 * 1.e0
- bfun(6) = sval**5 * 1.e0
- bfun(7) = sval**6 * 1.e0
+ bfun(1) = 1.e0_8
+ bfun(2) = sval**1 * 1.e0_8
+ bfun(3) = sval**2 * 1.e0_8
+ bfun(4) = sval**3 * 1.e0_8
+ bfun(5) = sval**4 * 1.e0_8
+ bfun(6) = sval**5 * 1.e0_8
+ bfun(7) = sval**6 * 1.e0_8
end select
@@ -153,46 +153,46 @@
!------------------------------------ +1st-order basis !
select case (ndof)
case (+1)
- bfun(1) = 0.e0
+ bfun(1) = 0.e0_8
case (+2)
- bfun(1) = 0.e0
- bfun(2) = 1.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 1.e0_8
case (+3)
- bfun(1) = 0.e0
- bfun(2) = 1.e0
- bfun(3) = sval**1 * 2.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 1.e0_8
+ bfun(3) = sval**1 * 2.e0_8
case (+4)
- bfun(1) = 0.e0
- bfun(2) = 1.e0
- bfun(3) = sval**1 * 2.e0
- bfun(4) = sval**2 * 3.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 1.e0_8
+ bfun(3) = sval**1 * 2.e0_8
+ bfun(4) = sval**2 * 3.e0_8
case (+5)
- bfun(1) = 0.e0
- bfun(2) = 1.e0
- bfun(3) = sval**1 * 2.e0
- bfun(4) = sval**2 * 3.e0
- bfun(5) = sval**3 * 4.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 1.e0_8
+ bfun(3) = sval**1 * 2.e0_8
+ bfun(4) = sval**2 * 3.e0_8
+ bfun(5) = sval**3 * 4.e0_8
case (+6)
- bfun(1) = 0.e0
- bfun(2) = 1.e0
- bfun(3) = sval**1 * 2.e0
- bfun(4) = sval**2 * 3.e0
- bfun(5) = sval**3 * 4.e0
- bfun(6) = sval**4 * 5.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 1.e0_8
+ bfun(3) = sval**1 * 2.e0_8
+ bfun(4) = sval**2 * 3.e0_8
+ bfun(5) = sval**3 * 4.e0_8
+ bfun(6) = sval**4 * 5.e0_8
case (+7)
- bfun(1) = 0.e0
- bfun(2) = 1.e0
- bfun(3) = sval**1 * 2.e0
- bfun(4) = sval**2 * 3.e0
- bfun(5) = sval**3 * 4.e0
- bfun(6) = sval**4 * 5.e0
- bfun(7) = sval**5 * 6.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 1.e0_8
+ bfun(3) = sval**1 * 2.e0_8
+ bfun(4) = sval**2 * 3.e0_8
+ bfun(5) = sval**3 * 4.e0_8
+ bfun(6) = sval**4 * 5.e0_8
+ bfun(7) = sval**5 * 6.e0_8
end select
@@ -200,46 +200,46 @@
!------------------------------------ +2nd-order basis !
select case (ndof)
case (+1)
- bfun(1) = 0.e0
+ bfun(1) = 0.e0_8
case (+2)
- bfun(1) = 0.e0
- bfun(2) = 0.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 0.e0_8
case (+3)
- bfun(1) = 0.e0
- bfun(2) = 0.e0
- bfun(3) = 2.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 0.e0_8
+ bfun(3) = 2.e0_8
case (+4)
- bfun(1) = 0.e0
- bfun(2) = 0.e0
- bfun(3) = 2.e0
- bfun(4) = sval**1 * 6.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 0.e0_8
+ bfun(3) = 2.e0_8
+ bfun(4) = sval**1 * 6.e0_8
case (+5)
- bfun(1) = 0.e0
- bfun(2) = 0.e0
- bfun(3) = 2.e0
- bfun(4) = sval**1 * 6.e0
- bfun(5) = sval**2 *12.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 0.e0_8
+ bfun(3) = 2.e0_8
+ bfun(4) = sval**1 * 6.e0_8
+ bfun(5) = sval**2 *12.e0_8
case (+6)
- bfun(1) = 0.e0
- bfun(2) = 0.e0
- bfun(3) = 2.e0
- bfun(4) = sval**1 * 6.e0
- bfun(5) = sval**2 *12.e0
- bfun(6) = sval**3 *20.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 0.e0_8
+ bfun(3) = 2.e0_8
+ bfun(4) = sval**1 * 6.e0_8
+ bfun(5) = sval**2 *12.e0_8
+ bfun(6) = sval**3 *20.e0_8
case (+7)
- bfun(1) = 0.e0
- bfun(2) = 0.e0
- bfun(3) = 2.e0
- bfun(4) = sval**1 * 6.e0
- bfun(5) = sval**2 *12.e0
- bfun(6) = sval**3 *20.e0
- bfun(7) = sval**4 *30.e0
+ bfun(1) = 0.e0_8
+ bfun(2) = 0.e0_8
+ bfun(3) = 2.e0_8
+ bfun(4) = sval**1 * 6.e0_8
+ bfun(5) = sval**2 *12.e0_8
+ bfun(6) = sval**3 *20.e0_8
+ bfun(7) = sval**4 *30.e0_8
end select
diff --git a/ext/PPR/src/ffsl1d.h90 b/ext/PPR/src/ffsl1d.h90
index c2a3d9fa438ce31622a1c4cf17391e5f4ee350b7..45b140740af61bcfc6180e43774b324d121e93b0 100644
--- a/ext/PPR/src/ffsl1d.h90
+++ b/ext/PPR/src/ffsl1d.h90
@@ -250,7 +250,7 @@
do ipos = +2 , npos - 1
- if (uvel(ipos) .gt. +0.e0) then
+ if (uvel(ipos) .gt. +0.e0_8) then
!----------- integrate profile over upwind cell IPOS-1 !
@@ -264,8 +264,8 @@
& * tDEL / SPAC( +1)
end if
- ss11 = +1.e0 - 2.e0 * uCFL
- ss22 = +1.e0
+ ss11 = +1.e0_8 - 2.e0_8 * uCFL
+ ss22 = +1.e0_8
call bfun1d(-1,mdof,ss11,vv11)
call bfun1d(-1,mdof,ss22,vv22)
@@ -285,7 +285,7 @@
end do
else &
- & if (uvel(ipos) .lt. -0.e0) then
+ & if (uvel(ipos) .lt. -0.e0_8) then
!----------- integrate profile over upwind cell IPOS+0 !
@@ -299,8 +299,8 @@
& * tDEL / SPAC( +1)
end if
- ss11 = -1.e0 - 2.e0 * uCFL
- ss22 = -1.e0
+ ss11 = -1.e0_8 - 2.e0_8 * uCFL
+ ss22 = -1.e0_8
call bfun1d(-1,mdof,ss11,vv11)
call bfun1d(-1,mdof,ss22,vv22)
diff --git a/ext/PPR/src/inv.h90 b/ext/PPR/src/inv.h90
index daceabb9a808cc6348e29c44604b21daedda93e4..491a1192322ec33479ce35ec93e9f73aa8f7f2f8 100644
--- a/ext/PPR/src/inv.h90
+++ b/ext/PPR/src/inv.h90
@@ -139,7 +139,7 @@
!-------------------------------- C = C + scal * A * B !
- if (scal .eq. +1.e0) then
+ if (scal .eq. +1.e0_8) then
cmat(1,1) = cmat(1,1) &
+ ( amat(1,1) * bmat(1,1) &
@@ -156,7 +156,7 @@
+ amat(2,2) * bmat(2,2) )
else &
- if (scal .eq. -1.e0) then
+ if (scal .eq. -1.e0_8) then
cmat(1,1) = cmat(1,1) &
- ( amat(1,1) * bmat(1,1) &
@@ -210,7 +210,7 @@
!-------------------------------- C = C + scal * A * B !
- if (scal .eq. +1.e0) then
+ if (scal .eq. +1.e0_8) then
cmat(1,1) = cmat(1,1) &
+ ( amat(1,1) * bmat(1,1) &
@@ -252,7 +252,7 @@
+ amat(3,3) * bmat(3,3) )
else &
- if (scal .eq. -1.e0) then
+ if (scal .eq. -1.e0_8) then
cmat(1,1) = cmat(1,1) &
- ( amat(1,1) * bmat(1,1) &
@@ -499,23 +499,23 @@
!---------------------------------------- L = C * A^-1 !
- lmat(1,1) = +0.e0
- lmat(1,2) = +0.e0
- lmat(2,1) = +0.e0
- lmat(2,2) = +0.e0
+ lmat(1,1) = +0.e0_8
+ lmat(1,2) = +0.e0_8
+ lmat(2,1) = +0.e0_8
+ lmat(2,2) = +0.e0_8
call mul_2x2(amat(3,1),adim,ainv,LDIM, &
- +1.e0,lmat,LDIM)
+ +1.e0_8,lmat,LDIM)
!---------------------------------------- U = A^-1 * B !
- umat(1,1) = +0.e0
- umat(1,2) = +0.e0
- umat(2,1) = +0.e0
- umat(2,2) = +0.e0
+ umat(1,1) = +0.e0_8
+ umat(1,2) = +0.e0_8
+ umat(2,1) = +0.e0_8
+ umat(2,2) = +0.e0_8
call mul_2x2(ainv,LDIM,amat(1,3),adim, &
- +1.e0,umat,LDIM)
+ +1.e0_8,umat,LDIM)
!-------------------------------- S = D - C * A^-1 * B !
@@ -525,7 +525,7 @@
smat(2,2) = amat(4,4)
call mul_2x2(lmat,LDIM,amat(1,3),adim, &
- -1.e0/adet,smat,LDIM)
+ -1.e0_8/adet,smat,LDIM)
call inv_2x2(smat,LDIM,sinv,LDIM,sdet)
@@ -640,33 +640,33 @@
!---------------------------------------- L = C * A^-1 !
- lmat(1,1) = +0.e0
- lmat(1,2) = +0.e0
- lmat(1,3) = +0.e0
- lmat(2,1) = +0.e0
- lmat(2,2) = +0.e0
- lmat(2,3) = +0.e0
- lmat(3,1) = +0.e0
- lmat(3,2) = +0.e0
- lmat(3,3) = +0.e0
+ lmat(1,1) = +0.e0_8
+ lmat(1,2) = +0.e0_8
+ lmat(1,3) = +0.e0_8
+ lmat(2,1) = +0.e0_8
+ lmat(2,2) = +0.e0_8
+ lmat(2,3) = +0.e0_8
+ lmat(3,1) = +0.e0_8
+ lmat(3,2) = +0.e0_8
+ lmat(3,3) = +0.e0_8
call mul_3x3(amat(4,1),adim,ainv,LDIM, &
- +1.e0,lmat,LDIM)
+ +1.e0_8,lmat,LDIM)
!---------------------------------------- U = A^-1 * B !
- umat(1,1) = +0.e0
- umat(1,2) = +0.e0
- umat(1,3) = +0.e0
- umat(2,1) = +0.e0
- umat(2,2) = +0.e0
- umat(2,3) = +0.e0
- umat(3,1) = +0.e0
- umat(3,2) = +0.e0
- umat(3,3) = +0.e0
+ umat(1,1) = +0.e0_8
+ umat(1,2) = +0.e0_8
+ umat(1,3) = +0.e0_8
+ umat(2,1) = +0.e0_8
+ umat(2,2) = +0.e0_8
+ umat(2,3) = +0.e0_8
+ umat(3,1) = +0.e0_8
+ umat(3,2) = +0.e0_8
+ umat(3,3) = +0.e0_8
call mul_3x3(ainv,LDIM,amat(1,4),adim, &
- +1.e0,umat,LDIM)
+ +1.e0_8,umat,LDIM)
!-------------------------------- S = D - C * A^-1 * B !
@@ -681,7 +681,7 @@
smat(3,3) = amat(6,6)
call mul_3x3(lmat,LDIM,amat(1,4),adim, &
- -1.e0/adet,smat,LDIM)
+ -1.e0_8/adet,smat,LDIM)
call inv_3x3(smat,LDIM,sinv,LDIM,sdet)
diff --git a/ext/PPR/src/oscl1d.h90 b/ext/PPR/src/oscl1d.h90
index 30f50982f55c0bb862f73334e0219a6102f4621a..178365709b2fe1245ff5ee2055fff8f3dd79802b 100644
--- a/ext/PPR/src/oscl1d.h90
+++ b/ext/PPR/src/oscl1d.h90
@@ -66,8 +66,8 @@
!------------------------------- at least 3 grid-cells !
do ipos = +1, npos-1
do ivar = +1, nvar-0
- oscl(1,ivar,ipos) = +0.e0
- oscl(2,ivar,ipos) = +0.e0
+ oscl(1,ivar,ipos) = +0.e0_8
+ oscl(2,ivar,ipos) = +0.e0_8
end do
end do
end if
@@ -143,25 +143,25 @@
hhlc = hhll + hhcc
hhmm = hhll + hhcc + hhrr
- cmat(1,1) = -(hhcc+2.e0*hhrr)/(hhlc*hhmm)
+ cmat(1,1) = -(hhcc+2.e0_8*hhrr)/(hhlc*hhmm)
cmat(1,2) = -(hhll-hhrr)* &
- & (3.e0*hhcc+2.e0*(hhll+hhrr))/&
+ & (3.e0_8*hhcc+2.e0_8*(hhll+hhrr))/&
& (hhlc*hhrc*hhmm)
- cmat(1,3) = +(hhcc+2.e0*hhll)/(hhrc*hhmm)
+ cmat(1,3) = +(hhcc+2.e0_8*hhll)/(hhrc*hhmm)
- cmat(2,1) = +3.e0/(hhlc*hhmm)
- cmat(2,2) = -3.e0*(2.e0*hhcc+hhll+hhrr)/&
+ cmat(2,1) = +3.e0_8/(hhlc*hhmm)
+ cmat(2,2) = -3.e0_8*(2.e0_8*hhcc+hhll+hhrr)/&
& (hhlc*hhrc*hhmm)
- cmat(2,3) = +3.e0/(hhrc*hhmm)
+ cmat(2,3) = +3.e0_8/(hhrc*hhmm)
do ivar = 1, nvar
- oscl(1,ivar,ipos) = +1.e0 * ( &
+ oscl(1,ivar,ipos) = +1.e0_8 * ( &
& + cmat(1,1)*fdat(1,ivar,ipos-1) &
& + cmat(1,2)*fdat(1,ivar,ipos+0) &
& + cmat(1,3)*fdat(1,ivar,ipos+1) )
- oscl(2,ivar,ipos) = +2.e0 * ( &
+ oscl(2,ivar,ipos) = +2.e0_8 * ( &
& + cmat(2,1)*fdat(1,ivar,ipos-1) &
& + cmat(2,2)*fdat(1,ivar,ipos+0) &
& + cmat(2,3)*fdat(1,ivar,ipos+1) )
@@ -176,8 +176,8 @@
hhcc = max(delx(head+1),dmin)
hhrr = max(delx(head+2),dmin)
- cmat(1,1) = -2.e0 / (hhll+hhcc)
- cmat(1,2) = +2.e0 / (hhll+hhcc)
+ cmat(1,1) = -2.e0_8 / (hhll+hhcc)
+ cmat(1,2) = +2.e0_8 / (hhll+hhcc)
do ivar = 1, nvar
@@ -185,7 +185,7 @@
& + cmat(1,1)*fdat(1,ivar,head+0) &
& + cmat(1,2)*fdat(1,ivar,head+1)
- oscl(2,ivar,head) = +0.e0
+ oscl(2,ivar,head) = +0.e0_8
end do
@@ -195,8 +195,8 @@
hhcc = max(delx(tail-1),dmin)
hhrr = max(delx(tail-0),dmin)
- cmat(1,2) = -2.e0 / (hhrr+hhcc)
- cmat(1,3) = +2.e0 / (hhrr+hhcc)
+ cmat(1,2) = -2.e0_8 / (hhrr+hhcc)
+ cmat(1,3) = +2.e0_8 / (hhrr+hhcc)
do ivar = 1, nvar
@@ -204,7 +204,7 @@
& + cmat(1,2)*fdat(1,ivar,tail-1) &
& + cmat(1,3)*fdat(1,ivar,tail+0)
- oscl(2,ivar,tail) = +0.e0
+ oscl(2,ivar,tail) = +0.e0_8
end do
@@ -268,7 +268,7 @@
& + .50d+0 * fdat(1,ivar,head+1) &
& - .50d+0 * fdat(1,ivar,head+0)
- oscl(2,ivar,head) = +0.e0
+ oscl(2,ivar,head) = +0.e0_8
end do
@@ -280,7 +280,7 @@
& + .50d+0 * fdat(1,ivar,tail+0) &
& - .50d+0 * fdat(1,ivar,tail-1)
- oscl(2,ivar,tail) = +0.e0
+ oscl(2,ivar,tail) = +0.e0_8
end do
diff --git a/ext/PPR/src/p1e.h90 b/ext/PPR/src/p1e.h90
index 39b639357f49aceba928005a09220803a9a00830..467c3507f633ecdcadbe178ac532a353a7fad249 100644
--- a/ext/PPR/src/p1e.h90
+++ b/ext/PPR/src/p1e.h90
@@ -83,10 +83,10 @@
do ivar = 1,nvar
edge(ivar,1) = fdat(1,ivar,1)
- dfdx(ivar,1) = 0.e0
+ dfdx(ivar,1) = 0.e0_8
edge(ivar,2) = fdat(1,ivar,1)
- dfdx(ivar,2) = 0.e0
+ dfdx(ivar,2) = 0.e0_8
end do
end if
@@ -103,7 +103,7 @@
!--------------- reconstruction: constant grid-spacing !
- dd10 = delx(+1) * 2.e0
+ dd10 = delx(+1) * 2.e0_8
do ivar = +1, nvar
@@ -175,8 +175,8 @@
edge(ivar,tail+1) = &
& fdat(+1,ivar,tail+0)
- dfdx(ivar,head-1) = 0.e0
- dfdx(ivar,tail+1) = 0.e0
+ dfdx(ivar,head-1) = 0.e0_8
+ dfdx(ivar,tail+1) = 0.e0_8
end do
diff --git a/ext/PPR/src/p3e.h90 b/ext/PPR/src/p3e.h90
index 6becd63358c1a4c9d1538386b65f59e91d164aea..a651873433e6cfbbc90e9242501c086b4af7bb57 100644
--- a/ext/PPR/src/p3e.h90
+++ b/ext/PPR/src/p3e.h90
@@ -120,24 +120,24 @@
do ivar = 1, nvar
edge(ivar,ipos) = ( &
- & - 1.e0 * &
+ & - 1.e0_8 * &
& fdat(1,ivar,ipos-2) &
- & + 7.e0 * &
+ & + 7.e0_8 * &
& fdat(1,ivar,ipos-1) &
- & + 7.e0 * &
+ & + 7.e0_8 * &
& fdat(1,ivar,ipos+0) &
- & - 1.e0 * &
- & fdat(1,ivar,ipos+1) ) / 12.e0
+ & - 1.e0_8 * &
+ & fdat(1,ivar,ipos+1) ) / 12.e0_8
dfdx(ivar,ipos) = ( &
- & + 1.e0 * &
+ & + 1.e0_8 * &
& fdat(1,ivar,ipos-2) &
- & - 15.e0 * &
+ & - 15.e0_8 * &
& fdat(1,ivar,ipos-1) &
- & + 15.e0 * &
+ & + 15.e0_8 * &
& fdat(1,ivar,ipos+0) &
- & - 1.e0 * &
- & fdat(1,ivar,ipos+1) ) / 12.e0
+ & - 1.e0_8 * &
+ & fdat(1,ivar,ipos+1) ) / 12.e0_8
dfdx(ivar,ipos) = &
& dfdx(ivar,ipos) / delx(+1)
@@ -165,7 +165,7 @@
xmap(-2) = -( delh(-2) &
& + delh(-1) ) / xhat
xmap(-1) = - delh(-1) / xhat
- xmap(+0) = + 0.e0
+ xmap(+0) = + 0.e0_8
xmap(+1) = + delh(+0) / xhat
xmap(+2) = +( delh(+0) &
& + delh(+1) ) / xhat
diff --git a/ext/PPR/src/p5e.h90 b/ext/PPR/src/p5e.h90
index fe4b446f7c1beec2e8a120c8902510a805d3c421..1e3f65fe93b366ec1d75b1948232e19474bc66a0 100644
--- a/ext/PPR/src/p5e.h90
+++ b/ext/PPR/src/p5e.h90
@@ -120,31 +120,31 @@
do ivar = 1, nvar
edge(ivar,ipos) = &
- & + ( 1.e0 / 60.e0) * &
+ & + ( 1.e0_8 / 60.e0_8) * &
& fdat(1,ivar,ipos-3) &
- & - ( 8.e0 / 60.e0) * &
+ & - ( 8.e0_8 / 60.e0_8) * &
& fdat(1,ivar,ipos-2) &
- & + (37.e0 / 60.e0) * &
+ & + (37.e0_8 / 60.e0_8) * &
& fdat(1,ivar,ipos-1) &
- & + (37.e0 / 60.e0) * &
+ & + (37.e0_8 / 60.e0_8) * &
& fdat(1,ivar,ipos+0) &
- & - ( 8.e0 / 60.e0) * &
+ & - ( 8.e0_8 / 60.e0_8) * &
& fdat(1,ivar,ipos+1) &
- & + ( 1.e0 / 60.e0) * &
+ & + ( 1.e0_8 / 60.e0_8) * &
& fdat(1,ivar,ipos+2)
dfdx(ivar,ipos) = &
- & - ( 1.e0 / 90.e0) * &
+ & - ( 1.e0_8 / 90.e0_8) * &
& fdat(1,ivar,ipos-3) &
- & + ( 5.e0 / 36.e0) * &
+ & + ( 5.e0_8 / 36.e0_8) * &
& fdat(1,ivar,ipos-2) &
- & - (49.e0 / 36.e0) * &
+ & - (49.e0_8 / 36.e0_8) * &
& fdat(1,ivar,ipos-1) &
- & + (49.e0 / 36.e0) * &
+ & + (49.e0_8 / 36.e0_8) * &
& fdat(1,ivar,ipos+0) &
- & - ( 5.e0 / 36.e0) * &
+ & - ( 5.e0_8 / 36.e0_8) * &
& fdat(1,ivar,ipos+1) &
- & + ( 1.e0 / 90.e0) * &
+ & + ( 1.e0_8 / 90.e0_8) * &
& fdat(1,ivar,ipos+2)
dfdx(ivar,ipos) = &
@@ -184,7 +184,7 @@
xmap(-2) = -( delh(-2) &
& + delh(-1) ) / xhat
xmap(-1) = - delh(-1) / xhat
- xmap(+0) = + 0.e0
+ xmap(+0) = + 0.e0_8
xmap(+1) = + delh(+0) / xhat
xmap(+2) = +( delh(+0) &
& + delh(+1) ) / xhat
diff --git a/ext/PPR/src/pbc.h90 b/ext/PPR/src/pbc.h90
index 1d09e3a9bc810bec8824ca176c579f7117e94339..13debd46f9e7adc0162e57b36b5c8e920f473370 100644
--- a/ext/PPR/src/pbc.h90
+++ b/ext/PPR/src/pbc.h90
@@ -244,8 +244,8 @@
xmap(-1) =-(delh(-1) + &
& delh(+0)*0.5d0)/xhat
- xmap(+0) = -1.e0
- xmap(+1) = +1.e0
+ xmap(+0) = -1.e0_8
+ xmap(+1) = +1.e0_8
xmap(+2) = (delh(+1) + &
& delh(+0)*0.5d0)/xhat
@@ -453,7 +453,7 @@
if (bcon(ivar)%bcopt.eq.bopt) then
eval(-1) = &
- & fdat(1,ivar,head-1) * 1.e0
+ & fdat(1,ivar,head-1) * 1.e0_8
eval(+0) = &
& fdat(1,ivar,head-1) * .5d0 + &
& fdat(1,ivar,head+0) * .5d0
@@ -575,8 +575,8 @@
xmap(-1) =-(delh(-1) + &
& delh(+0)*0.5d0)/xhat
- xmap(+0) = -1.e0
- xmap(+1) = +1.e0
+ xmap(+0) = -1.e0_8
+ xmap(+1) = +1.e0_8
xmap(+2) = (delh(+1) + &
& delh(+0)*0.5d0)/xhat
@@ -790,7 +790,7 @@
& fdat(1,ivar,tail+0) * .5d0 + &
& fdat(1,ivar,tail+1) * .5d0
eval(+2) = &
- & fdat(1,ivar,tail+1) * 1.e0
+ & fdat(1,ivar,tail+1) * 1.e0_8
gval(+0) = &
& fdat(1,ivar,tail+0) * .5d0 - &
diff --git a/ext/PPR/src/ppm.h90 b/ext/PPR/src/ppm.h90
index f8158111cfd1c1e9b6a285e4ed357fe433e51124..923c80f956f0e3f296a11efca8380d84397127eb 100644
--- a/ext/PPR/src/ppm.h90
+++ b/ext/PPR/src/ppm.h90
@@ -94,8 +94,8 @@
do ivar = +1, nvar
fhat(1,ivar,+1) = &
& fdat(1,ivar,+1)
- fhat(2,ivar,+1) = 0.e0
- fhat(3,ivar,+1) = 0.e0
+ fhat(2,ivar,+1) = 0.e0_8
+ fhat(3,ivar,+1) = 0.e0_8
end do
end if
@@ -275,8 +275,8 @@
mono = +1
lhat(1) = ff00
- lhat(2) = 0.e0
- lhat(3) = 0.e0
+ lhat(2) = 0.e0_8
+ lhat(3) = 0.e0_8
return
diff --git a/ext/PPR/src/pqm.h90 b/ext/PPR/src/pqm.h90
index c834c3daa5886871c672841bd99239ce62c09ce2..19db6503c361f0e99950fb726fae12fa237cec60 100644
--- a/ext/PPR/src/pqm.h90
+++ b/ext/PPR/src/pqm.h90
@@ -99,10 +99,10 @@
!----- default to reduced order if insufficient points !
do ivar = +1, nvar
fhat(1,ivar,+1) = fdat(1,ivar,+1)
- fhat(2,ivar,+1) = 0.e0
- fhat(3,ivar,+1) = 0.e0
- fhat(4,ivar,+1) = 0.e0
- fhat(5,ivar,+1) = 0.e0
+ fhat(2,ivar,+1) = 0.e0_8
+ fhat(3,ivar,+1) = 0.e0_8
+ fhat(4,ivar,+1) = 0.e0_8
+ fhat(5,ivar,+1) = 0.e0_8
end do
end if
@@ -313,10 +313,10 @@
mono = +1
lhat(1) = ff00
- lhat(2) = 0.e0
- lhat(3) = 0.e0
- lhat(4) = 0.e0
- lhat(5) = 0.e0
+ lhat(2) = 0.e0_8
+ lhat(3) = 0.e0_8
+ lhat(4) = 0.e0_8
+ lhat(5) = 0.e0_8
return
diff --git a/ext/PPR/src/rmap1d.h90 b/ext/PPR/src/rmap1d.h90
index fa3486e8b022199df66f9d71fb53aa52e3515c41..dc10cd74ee1530d390157a00135a916d2076ef52 100644
--- a/ext/PPR/src/rmap1d.h90
+++ b/ext/PPR/src/rmap1d.h90
@@ -244,8 +244,8 @@
!------------------------------------- initializations !
- vvlo(+1:+5) = 0.e0
- vvhi(+1:+5) = 0.e0
+ vvlo(+1:+5) = 0.e0_8
+ vvhi(+1:+5) = 0.e0_8
!------------- remap FDAT from XPOS to XNEW using FHAT !
@@ -285,7 +285,7 @@
do idof = +1,ndof
do ivar = +1,nvar
- fnew(idof,ivar,kpos) = 0.e0
+ fnew(idof,ivar,kpos) = 0.e0_8
end do
end do
@@ -374,7 +374,7 @@
! wrt. column-wise conservation. Use KBN approach to
! account for FP roundoff.
- sdat = 0.e0; serr = 0.e0
+ sdat = 0.e0_8; serr = 0.e0_8
do ipos = +1, npos-1
do ivar = +1, nvar-0
@@ -407,7 +407,7 @@
sdat = sdat + serr
- snew = 0.e0; serr = 0.e0
+ snew = 0.e0_8; serr = 0.e0_8
do ipos = +1, nnew-1
do ivar = +1, nvar-0
@@ -450,7 +450,7 @@
do ivar = +1, nvar-0
- if (serr(ivar) .gt. 0.e0) then
+ if (serr(ivar) .gt. 0.e0_8) then
vmin = kmin(ivar)
@@ -459,7 +459,7 @@
& serr(ivar)/(xnew(vmin+1)-xnew(vmin+0))
else &
- & if (serr(ivar) .lt. 0.e0) then
+ & if (serr(ivar) .lt. 0.e0_8) then
vmax = kmax(ivar)
diff --git a/ext/PPR/src/util1d.h90 b/ext/PPR/src/util1d.h90
index 6b4353ba76ee5c76aaf0541e8f191770ba0ecfd1..10cd4c58c573e501145b647580300ee00ff85160 100644
--- a/ext/PPR/src/util1d.h90
+++ b/ext/PPR/src/util1d.h90
@@ -110,7 +110,7 @@
call random_number (rand)
- rand = 2.e0 * (rand-.5d0)
+ rand = 2.e0_8 * (rand-.5d0)
move = rand * move
diff --git a/ext/PPR/src/weno1d.h90 b/ext/PPR/src/weno1d.h90
index d74c76faa23d142eba27ced7929cb5dc4e017c46..e65272ee7e9fd89709fc53cf1d0cb05fa94e1ca7 100644
--- a/ext/PPR/src/weno1d.h90
+++ b/ext/PPR/src/weno1d.h90
@@ -111,7 +111,7 @@
wsum = wval(1) + wval(2) + ZERO
wval(1) = wval(1) / wsum
! wval(2) = wval(2) / wsum
- wval(2) =-wval(1) + 1.e0 ! wval(2)/wsum but robust !
+ wval(2) =-wval(1) + 1.e0_8 ! wval(2)/wsum but robust !
return
@@ -168,15 +168,15 @@
if (ipos-halo.lt.head) then
- omax = 1.e0
- omin = 0.e0 ; return
+ omax = 1.e0_8
+ omin = 0.e0_8 ; return
end if
if (ipos+halo.gt.tail) then
- omax = 1.e0
- omin = 0.e0 ; return
+ omax = 1.e0_8
+ omin = 0.e0_8 ; return
end if
@@ -203,7 +203,7 @@
!---------------------------------------- "lower" part !
- delh = 0.e0
+ delh = 0.e0_8
do hpos = ipos-1, imin, -1
@@ -235,7 +235,7 @@
!---------------------------------------- "upper" part !
- delh = 0.e0
+ delh = 0.e0_8
do hpos = ipos+1, imax, +1
@@ -319,15 +319,15 @@
if (ipos-halo.lt.head) then
- omax = 1.e0
- omin = 0.e0 ; return
+ omax = 1.e0_8
+ omin = 0.e0_8 ; return
end if
if (ipos+halo.gt.tail) then
- omax = 1.e0
- omin = 0.e0 ; return
+ omax = 1.e0_8
+ omin = 0.e0_8 ; return
end if
@@ -343,21 +343,21 @@
dfx1 = oscl(1,ivar,ipos)
dfx2 = oscl(2,ivar,ipos)
- oval = (2.e0**1*dfx1)**2 &
- & + (2.e0**2*dfx2)**2
+ oval = (2.e0_8**1*dfx1)**2 &
+ & + (2.e0_8**2*dfx2)**2
omin = oval
omax = oval
!---------------------------------------- "lower" part !
- delh = 0.e0
+ delh = 0.e0_8
do hpos = ipos-1, imin, -1
!------------------ calc. derivatives centred on IPOS. !
- delh = delh + 2.e0
+ delh = delh + 2.e0_8
dfx1 = oscl(1,ivar,hpos)
dfx2 = oscl(2,ivar,hpos)
@@ -366,8 +366,8 @@
!------------------ indicator: NORM(H^N * D^N/DX^N(F)) !
- oval = (2.e0**1*dfx1)**2 &
- & + (2.e0**2*dfx2)**2
+ oval = (2.e0_8**1*dfx1)**2 &
+ & + (2.e0_8**2*dfx2)**2
if (oval .lt. omin) then
omin = oval
@@ -380,13 +380,13 @@
!---------------------------------------- "upper" part !
- delh = 0.e0
+ delh = 0.e0_8
do hpos = ipos+1, imax, +1
!------------------ calc. derivatives centred on IPOS. !
- delh = delh - 2.e0
+ delh = delh - 2.e0_8
dfx1 = oscl(1,ivar,hpos)
dfx2 = oscl(2,ivar,hpos)
@@ -395,8 +395,8 @@
!------------------ indicator: NORM(H^N * D^N/DX^N(F)) !
- oval = (2.e0**1*dfx1)**2 &
- & + (2.e0**2*dfx2)**2
+ oval = (2.e0_8**1*dfx1)**2 &
+ & + (2.e0_8**2*dfx2)**2
if (oval .lt. omin) then
omin = oval
diff --git a/src/ICE/icedyn_adv_umx.F90 b/src/ICE/icedyn_adv_umx.F90
index 2aa6a5eb7983c14b9485aeee34b4f18d4d70b41c..755c002cb80be72346220f854dc0c99ce8d5d8ae 100644
--- a/src/ICE/icedyn_adv_umx.F90
+++ b/src/ICE/icedyn_adv_umx.F90
@@ -1165,7 +1165,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj ) :: zbup, zbdo
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ups, ztj_ups
!!----------------------------------------------------------------------
- zbig = 1.e+40_wp
+ zbig = 1.e+20_wp ! works ok with simple/double precison
! antidiffusive flux : high order minus low order
! --------------------------------------------------
diff --git a/src/ICE/icedyn_rhg_vp.F90 b/src/ICE/icedyn_rhg_vp.F90
index 3486e77cd349c835d0f38f6c419b60cbe2f05ea9..81e636c94cb2f2aa859823483b90591d3b938bda 100644
--- a/src/ICE/icedyn_rhg_vp.F90
+++ b/src/ICE/icedyn_rhg_vp.F90
@@ -1124,7 +1124,7 @@ CONTAINS
END_2D
-!!$ CALL lbc_lnk( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. )
+!!$ CALL lbc_lnk( 'icedyn_rhg_vp', zs1, 'T', 1._wp, zs2, 'T', 1._wp, zs12, 'T', 1._wp )
ENDIF
@@ -1141,7 +1141,7 @@ CONTAINS
END_2D
- CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. )
+ CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1._wp )
ENDIF
@@ -1178,8 +1178,8 @@ CONTAINS
END_2D
!
- CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, &
-! & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
+ CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1._wp, ztauy_oi, 'V', -1._wp, ztaux_ai, 'U', -1._wp, ztauy_ai, 'V', -1._wp ) !, &
+! & ztaux_bi, 'U', -1._wp, ztauy_bi, 'V', -1._wp )
!
CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
@@ -1209,7 +1209,7 @@ CONTAINS
zsig_II(ji,jj) = 0.5_wp * SQRT ( zs2(ji,jj) * zs2(ji,jj) + 4. * zs12(ji,jj) * zs12(ji,jj) )
END_2D
-!!$ CALL lbc_lnk( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.)
+!!$ CALL lbc_lnk( 'icedyn_rhg_vp', zsig_I, 'T', 1._wp, zsig_II, 'T', 1._wp)
IF( iom_use('normstr') ) CALL iom_put( 'normstr' , zsig_I(:,:) * zmsk00(:,:) ) ! Normal stress
IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
@@ -1250,7 +1250,7 @@ CONTAINS
END_2D
!
- ! CALL lbc_lnk( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
+ ! CALL lbc_lnk( 'icedyn_rhg_vp', zsig1_p, 'T', 1._wp, zsig2_p, 'T', 1._wp)
!
CALL iom_put( 'sig1_pnorm' , zsig1_p )
CALL iom_put( 'sig2_pnorm' , zsig2_p )
@@ -1274,8 +1274,8 @@ CONTAINS
& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
END_2D
!
- CALL lbc_lnk( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., &
- & zCorU, 'U', -1., zCorV, 'V', -1. )
+ CALL lbc_lnk( 'icedyn_rhg_vp', zspgU, 'U', -1._wp, zspgV, 'V', -1._wp, &
+ & zCorU, 'U', -1._wp, zCorV, 'V', -1._wp )
!
CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x)
CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y)
@@ -1305,7 +1305,7 @@ CONTAINS
END_2D
- CALL lbc_lnk( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. )
+ CALL lbc_lnk( 'icedyn_rhg_vp', zfU, 'U', -1._wp, zfV, 'V', -1._wp )
CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x)
CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y)
@@ -1335,9 +1335,9 @@ CONTAINS
END_2D
- CALL lbc_lnk( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
- & zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
- & zdiag_xatrp , 'U', -1., zdiag_yatrp , 'V', -1. )
+ CALL lbc_lnk( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1._wp, zdiag_ymtrp_ice, 'V', -1._wp, &
+ & zdiag_xmtrp_snw, 'U', -1._wp, zdiag_ymtrp_snw, 'V', -1._wp, &
+ & zdiag_xatrp , 'U', -1._wp, zdiag_yatrp , 'V', -1._wp )
CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s)
CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport
@@ -1628,7 +1628,7 @@ CONTAINS
pvel_res(ji,jj) = 0.25_wp * ( ( zu_res(ji-1,jj) + zu_res(ji,jj) ) + ( zv_res(ji,jj-1) + zv_res(ji,jj) ) )
END_2D
- CALL lbc_lnk( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1. )
+ CALL lbc_lnk( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1._wp )
ELSE
diff --git a/src/ICE/icethd.F90 b/src/ICE/icethd.F90
index ddc4ddf2b6fdce2307bbf7f398c74e17f34c0b17..8408ab9891c9e055cb94a327629ef66a588e2461 100644
--- a/src/ICE/icethd.F90
+++ b/src/ICE/icethd.F90
@@ -345,11 +345,22 @@ CONTAINS
CALL tab_2d_1d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
!
! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
+ ! Here we make sure that we don't divide by very small, but physically
+ ! meaningless, products of sea ice thicknesses/snow depths and sea ice
+ ! concentration
DO jk = 1, nlay_i
- WHERE( h_i_1d(1:npti)>0._wp ) e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i
+ WHERE( (h_i_1d(1:npti) * a_i_1d(1:npti)) > epsi20 )
+ e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i
+ ELSEWHERE
+ e_i_1d(1:npti,jk) = 0._wp
+ ENDWHERE
END DO
DO jk = 1, nlay_s
- WHERE( h_s_1d(1:npti)>0._wp ) e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s
+ WHERE( (h_s_1d(1:npti) * a_i_1d(1:npti)) > epsi20 )
+ e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s
+ ELSEWHERE
+ e_s_1d(1:npti,jk) = 0._wp
+ ENDWHERE
END DO
!
! !---------------------!
diff --git a/src/NST/vremap.F90 b/src/NST/vremap.F90
index 6e3853bff00e1117653109e5c839ae2d31aab04e..a16b1f83d21dd9e6666c5ebe9f6e0faf26b9c0fb 100644
--- a/src/NST/vremap.F90
+++ b/src/NST/vremap.F90
@@ -292,8 +292,8 @@ CONTAINS
!
INTEGER, PARAMETER :: ndof = 1
INTEGER :: jk, jn
- REAL(wp) :: zwin(kjpk_in+1) , ztin(ndof, kn_var, kjpk_in)
- REAL(wp) :: zwout(kjpk_out+1), ztout(ndof, kn_var, kjpk_out)
+ REAL(dp) :: zwin(kjpk_in+1) , ztin(ndof, kn_var, kjpk_in) ! rmap1d uses dp
+ REAL(dp) :: zwout(kjpk_out+1), ztout(ndof, kn_var, kjpk_out) ! rmap1d uses dp
TYPE(rmap_work) :: work
TYPE(rmap_opts) :: opts
TYPE(rcon_ends) :: bc_l(kn_var)
diff --git a/src/OCE/DOM/domzgr.F90 b/src/OCE/DOM/domzgr.F90
index a151f2f46ab25b4a05bfdbe81b034e5a5d32f0b5..9029c88e79b3b0c632127c5457550a0b69de71df 100644
--- a/src/OCE/DOM/domzgr.F90
+++ b/src/OCE/DOM/domzgr.F90
@@ -293,7 +293,7 @@ CALL FLUSH(numout)
IF( .NOT. ( l_Jperio .OR. l_NFold ) ) THEN ! N closed:
zmsk(:,mj0(jpjglo-nn_hls,nn_hls):mj1(jpjglo-nn_hls,nn_hls) ) = 0._wp ! last line of inner global domain at 0
ENDIF
- CALL lbc_lnk( 'usrdef_zgr', zmsk, 'T', 1. ) ! set halos
+ CALL lbc_lnk( 'usrdef_zgr', zmsk, 'T', 1._wp ) ! set halos
k_top(:,:) = k_top(:,:) * NINT( zmsk(:,:) )
!
! Any closed seas (defined by closea_mask > 0 in domain_cfg file) to be filled
diff --git a/src/OCE/IOM/iom.F90 b/src/OCE/IOM/iom.F90
index 6244df8c386c83acc0df99ac11ce873ba98783d7..382e158778991789c23b4adb1d0ab732ad4e975e 100644
--- a/src/OCE/IOM/iom.F90
+++ b/src/OCE/IOM/iom.F90
@@ -363,6 +363,7 @@ CONTAINS
CHARACTER(LEN=nf90_max_name), ALLOCATABLE :: indimnames(:)
CHARACTER(LEN=nf90_max_name) :: dimname, varname
INTEGER :: iln
+ INTEGER :: iprec
CHARACTER(LEN=lc) :: fname
LOGICAL :: lmeta
!metadata in restart file for restart read with XIOS
@@ -385,7 +386,7 @@ CONTAINS
clinfo = ' iom_set_vars_active, file: '//TRIM(iom_file(idnum)%name)
iln = INDEX( iom_file(idnum)%name, '.nc' )
-!XIOS doee not need .nc
+!XIOS does not need .nc
IF(iln > 0) THEN
fname = iom_file(idnum)%name(1:iln-1)
ELSE
@@ -419,29 +420,33 @@ CONTAINS
ENDDO
IF(.NOT.lmeta) THEN
CALL xios_add_child(file_hdl, field_hdl, varname)
- mdims = ndims
- IF(ANY(dimids(1:ndims) == unlimitedDimId)) THEN
- mdims = mdims - 1
+ IF( xtype == NF90_FLOAT ) THEN ; iprec = sp
+ ELSEIF( xtype == NF90_DOUBLE ) THEN ; iprec = dp
+ ELSE
+ WRITE(ctmp1,*) 'iom_set_vars_active: variable ', TRIM(varname) ,' wrong type : ', xtype
+ CALL ctl_stop( 'iom_set_vars_active:', ctmp1 )
ENDIF
+
+ mdims = ndims - COUNT( (/ ANY(dimids(1:ndims) == unlimitedDimId) /) )
- IF(mdims == 3) THEN
- CALL xios_set_attr (field_hdl, enabled = .TRUE., name = varname, &
- domain_ref="grid_N", &
- axis_ref=iom_axis(indimlens(dimids(mdims))), &
- prec = 8, operation = "instant" )
+ IF( mdims == 3) THEN
+ CALL xios_set_attr (field_hdl, enabled = .TRUE., name = varname, &
+ domain_ref="grid_N", &
+ axis_ref=iom_axis(indimlens(dimids(mdims))), &
+ prec = iprec, operation = "instant" )
ELSEIF(mdims == 2) THEN
- CALL xios_set_attr (field_hdl, enabled = .TRUE., name = varname, &
- domain_ref="grid_N", prec = 8, &
- operation = "instant" )
+ CALL xios_set_attr (field_hdl, enabled = .TRUE., name = varname, &
+ domain_ref="grid_N", &
+ prec = iprec, operation = "instant" )
ELSEIF(mdims == 1) THEN
CALL xios_set_attr (field_hdl, enabled = .TRUE., name = varname, &
axis_ref=iom_axis(indimlens(dimids(mdims))), &
- prec = 8, operation = "instant" )
+ prec = iprec, operation = "instant" )
ELSEIF(mdims == 0) THEN
CALL xios_set_attr (field_hdl, enabled = .TRUE., name = varname, &
- scalar_ref = "grid_scalar", prec = 8, &
- operation = "instant" )
+ scalar_ref = "grid_scalar", &
+ prec = iprec, operation = "instant" )
ELSE
WRITE(ctmp1,*) 'iom_set_vars_active: variable ', TRIM(varname) ,' incorrect number of dimensions'
CALL ctl_stop( 'iom_set_vars_active:', ctmp1 )
@@ -534,40 +539,40 @@ CONTAINS
CALL ctl_stop(ctmp1, ctmp2)
ENDIF
- IF(PRESENT(rd3)) THEN
+ IF( PRESENT(rd3)) THEN
CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
domain_ref = "grid_N"//TRIM(clgsuf), &
axis_ref = iom_axis(size(rd3, 3)), &
- prec = 8, operation = "instant" )
+ prec = dp, operation = "instant" )
ELSEIF(PRESENT(rs3)) THEN
CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
domain_ref = "grid_N"//TRIM(clgsuf), &
axis_ref = iom_axis(size(rd3, 3)), &
- prec = 4, operation = "instant" )
+ prec = sp, operation = "instant" )
ELSEIF(PRESENT(rd2)) THEN
CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
- domain_ref = "grid_N"//TRIM(clgsuf), prec = 8, &
- operation = "instant" )
+ domain_ref = "grid_N"//TRIM(clgsuf), &
+ prec = dp, operation = "instant" )
ELSEIF(PRESENT(rs2)) THEN
CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
- domain_ref = "grid_N"//TRIM(clgsuf), prec = 4, &
- operation = "instant" )
+ domain_ref = "grid_N"//TRIM(clgsuf), &
+ prec = sp, operation = "instant" )
ELSEIF(PRESENT(rd1)) THEN
CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
axis_ref = iom_axis(size(rd1, 1)), &
- prec = 8, operation = "instant" )
+ prec = dp, operation = "instant" )
ELSEIF(PRESENT(rs1)) THEN
CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
axis_ref = iom_axis(size(rd1, 1)), &
- prec = 4, operation = "instant" )
+ prec = sp, operation = "instant" )
ELSEIF(PRESENT(rd0)) THEN
CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
- scalar_ref = "grid_scalar", prec = 8, &
- operation = "instant" )
+ scalar_ref = "grid_scalar", &
+ prec = dp, operation = "instant" )
ELSEIF(PRESENT(rs0)) THEN
CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
- scalar_ref = "grid_scalar", prec = 4, &
- operation = "instant" )
+ scalar_ref = "grid_scalar", &
+ prec = sp, operation = "instant" )
ENDIF
#endif
END SUBROUTINE iom_set_rstw_active
@@ -1077,19 +1082,13 @@ CONTAINS
INTEGER , INTENT(in ) :: kdom ! Type of domain to be read
CHARACTER(len=*), INTENT(in ) :: cdvar ! Name of the variable
REAL(sp) , INTENT( out), DIMENSION(:) :: pvar ! read field
- REAL(dp) , ALLOCATABLE , DIMENSION(:) :: ztmp_pvar ! tmp var to read field
INTEGER , INTENT(in ) , OPTIONAL :: ktime ! record number
INTEGER , INTENT(in ), DIMENSION(1), OPTIONAL :: kstart ! start axis position of the reading
INTEGER , INTENT(in ), DIMENSION(1), OPTIONAL :: kcount ! number of points in each axis
!
IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) THEN
- ALLOCATE(ztmp_pvar(size(pvar,1)))
- CALL iom_get_123d( kiomid, kdom , cdvar , pv_r1d=ztmp_pvar, &
- & ktime=ktime, kstart=kstart, kcount=kcount )
- pvar = ztmp_pvar
- DEALLOCATE(ztmp_pvar)
- END IF
+ IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvsp1d = pvar, &
+ & ktime = ktime, kstart = kstart, kcount = kcount )
ENDIF
END SUBROUTINE iom_g1d_sp
@@ -1104,8 +1103,8 @@ CONTAINS
INTEGER , INTENT(in ), DIMENSION(1), OPTIONAL :: kcount ! number of points in each axis
!
IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom , cdvar , pv_r1d=pvar, &
- & ktime=ktime, kstart=kstart, kcount=kcount)
+ IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvdp1d = pvar, &
+ & ktime = ktime, kstart = kstart, kcount = kcount)
ENDIF
END SUBROUTINE iom_g1d_dp
@@ -1114,23 +1113,17 @@ CONTAINS
INTEGER , INTENT(in ) :: kdom ! Type of domain to be read
CHARACTER(len=*), INTENT(in ) :: cdvar ! Name of the variable
REAL(sp) , INTENT( out), DIMENSION(:,:) :: pvar ! read field
- REAL(dp) , ALLOCATABLE , DIMENSION(:,:) :: ztmp_pvar ! tmp var to read field
INTEGER , INTENT(in ) , OPTIONAL :: ktime ! record number
CHARACTER(len=1), INTENT(in ) , OPTIONAL :: cd_type ! nature of grid-points (T, U, V, F, W)
- REAL(dp) , INTENT(in ) , OPTIONAL :: psgn ! -1.(1.): (not) change sign across the north fold
+ REAL(sp) , INTENT(in ) , OPTIONAL :: psgn ! -1.(1.): (not) change sign across the north fold
INTEGER , INTENT(in ) , OPTIONAL :: kfill ! value of kfillmode in lbc_lbk
INTEGER , INTENT(in ), DIMENSION(2), OPTIONAL :: kstart ! start axis position of the reading
INTEGER , INTENT(in ), DIMENSION(2), OPTIONAL :: kcount ! number of points in each axis
!
IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) THEN
- ALLOCATE(ztmp_pvar(size(pvar,1), size(pvar,2)))
- CALL iom_get_123d( kiomid, kdom, cdvar , pv_r2d = ztmp_pvar , ktime = ktime, &
- & cd_type = cd_type, psgn = psgn , kfill = kfill, &
- & kstart = kstart , kcount = kcount )
- pvar = ztmp_pvar
- DEALLOCATE(ztmp_pvar)
- ENDIF
+ IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvsp2d = pvar, &
+ & cd_type = cd_type, psgn_sp = psgn, kfill = kfill, &
+ & ktime = ktime, kstart = kstart, kcount = kcount )
ENDIF
END SUBROUTINE iom_g2d_sp
@@ -1147,9 +1140,9 @@ CONTAINS
INTEGER , INTENT(in ), DIMENSION(2), OPTIONAL :: kcount ! number of points in each axis
!
IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar , pv_r2d = pvar , ktime = ktime, &
- & cd_type = cd_type, psgn = psgn , kfill = kfill, &
- & kstart = kstart , kcount = kcount )
+ IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvdp2d = pvar, &
+ & cd_type = cd_type, psgn_dp = psgn, kfill = kfill, &
+ & ktime = ktime, kstart = kstart, kcount = kcount )
ENDIF
END SUBROUTINE iom_g2d_dp
@@ -1158,23 +1151,17 @@ CONTAINS
INTEGER , INTENT(in ) :: kdom ! Type of domain to be read
CHARACTER(len=*), INTENT(in ) :: cdvar ! Name of the variable
REAL(sp) , INTENT( out), DIMENSION(:,:,:) :: pvar ! read field
- REAL(dp) , ALLOCATABLE , DIMENSION(:,:,:) :: ztmp_pvar ! tmp var to read field
INTEGER , INTENT(in ) , OPTIONAL :: ktime ! record number
CHARACTER(len=1), INTENT(in ) , OPTIONAL :: cd_type ! nature of grid-points (T, U, V, F, W)
- REAL(dp) , INTENT(in ) , OPTIONAL :: psgn ! -1.(1.) : (not) change sign across the north fold
+ REAL(sp) , INTENT(in ) , OPTIONAL :: psgn ! -1.(1.) : (not) change sign across the north fold
INTEGER , INTENT(in ) , OPTIONAL :: kfill ! value of kfillmode in lbc_lbk
INTEGER , INTENT(in ), DIMENSION(3), OPTIONAL :: kstart ! start axis position of the reading
INTEGER , INTENT(in ), DIMENSION(3), OPTIONAL :: kcount ! number of points in each axis
!
IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) THEN
- ALLOCATE(ztmp_pvar(size(pvar,1), size(pvar,2), size(pvar,3)))
- CALL iom_get_123d( kiomid, kdom, cdvar , pv_r3d = ztmp_pvar , ktime = ktime, &
- & cd_type = cd_type, psgn = psgn , kfill = kfill, &
- & kstart = kstart , kcount = kcount )
- pvar = ztmp_pvar
- DEALLOCATE(ztmp_pvar)
- END IF
+ IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvsp3d = pvar, &
+ & cd_type = cd_type, psgn_sp = psgn, kfill = kfill, &
+ & ktime = ktime, kstart = kstart, kcount = kcount )
ENDIF
END SUBROUTINE iom_g3d_sp
@@ -1191,18 +1178,16 @@ CONTAINS
INTEGER , INTENT(in ), DIMENSION(3), OPTIONAL :: kcount ! number of points in each axis
!
IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) THEN
- CALL iom_get_123d( kiomid, kdom, cdvar , pv_r3d = pvar , ktime = ktime, &
- & cd_type = cd_type, psgn = psgn , kfill = kfill, &
- & kstart = kstart , kcount = kcount )
- END IF
+ IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvdp3d = pvar, &
+ & cd_type = cd_type, psgn_dp = psgn, kfill = kfill, &
+ & ktime = ktime, kstart = kstart, kcount = kcount )
ENDIF
END SUBROUTINE iom_g3d_dp
!!----------------------------------------------------------------------
- SUBROUTINE iom_get_123d( kiomid , kdom, cdvar, pv_r1d, pv_r2d, pv_r3d, ktime , &
- & cd_type, psgn, kfill, kstart, kcount )
+ SUBROUTINE iom_get_123d( kiomid , kdom, cdvar, pvsp1d, pvsp2d, pvsp3d, pvdp1d, pvdp2d, pvdp3d, &
+ & ktime , cd_type, psgn_sp, psgn_dp, kfill, kstart, kcount )
!!-----------------------------------------------------------------------
!! *** ROUTINE iom_get_123d ***
!!
@@ -1213,12 +1198,16 @@ CONTAINS
INTEGER , INTENT(in ) :: kiomid ! Identifier of the file
INTEGER , INTENT(in ) :: kdom ! Type of domain to be read
CHARACTER(len=*) , INTENT(in ) :: cdvar ! Name of the variable
- REAL(dp), DIMENSION(:) , INTENT( out), OPTIONAL :: pv_r1d ! read field (1D case)
- REAL(dp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pv_r2d ! read field (2D case)
- REAL(dp), DIMENSION(:,:,:) , INTENT( out), OPTIONAL :: pv_r3d ! read field (3D case)
+ REAL(sp), DIMENSION(:) , INTENT( out), OPTIONAL :: pvsp1d ! read field (1D case), single precision
+ REAL(sp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pvsp2d ! read field (2D case), single precision
+ REAL(sp), DIMENSION(:,:,:) , INTENT( out), OPTIONAL :: pvsp3d ! read field (3D case), single precision
+ REAL(dp), DIMENSION(:) , INTENT( out), OPTIONAL :: pvdp1d ! read field (1D case), double precision
+ REAL(dp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pvdp2d ! read field (2D case), double precision
+ REAL(dp), DIMENSION(:,:,:) , INTENT( out), OPTIONAL :: pvdp3d ! read field (3D case), double precision
INTEGER , INTENT(in ), OPTIONAL :: ktime ! record number
CHARACTER(len=1) , INTENT(in ), OPTIONAL :: cd_type ! nature of grid-points (T, U, V, F, W)
- REAL(dp) , INTENT(in ), OPTIONAL :: psgn ! -1.(1.) : (not) change sign across the north fold
+ REAL(sp) , INTENT(in ), OPTIONAL :: psgn_sp ! -1.(1.) : (not) change sign across the north fold
+ REAL(dp) , INTENT(in ), OPTIONAL :: psgn_dp ! -1.(1.) : (not) change sign across the north fold
INTEGER , INTENT(in ), OPTIONAL :: kfill ! value of kfillmode in lbc_lbk
INTEGER , DIMENSION(:) , INTENT(in ), OPTIONAL :: kstart ! start position of the reading in each axis
INTEGER , DIMENSION(:) , INTENT(in ), OPTIONAL :: kcount ! number of points to be read in each axis
@@ -1240,35 +1229,38 @@ CONTAINS
INTEGER, DIMENSION(jpmax_dims) :: icnt ! number of value to read along each axis
INTEGER, DIMENSION(jpmax_dims) :: idimsz ! size of the dimensions of the variable
INTEGER, DIMENSION(jpmax_dims) :: ishape ! size of the dimensions of the variable
- REAL(dp) :: zscf, zofs ! sacle_factor and add_offset
- REAL(wp) :: zsgn ! local value of psgn
+ REAL(sp) :: zscf_sp, zofs_sp ! sacle_factor and add_offset, single precision
+ REAL(dp) :: zscf_dp, zofs_dp ! sacle_factor and add_offset, double precision
+ REAL(sp) :: zsgn_sp ! local value of psgn, single precision
+ REAL(dp) :: zsgn_dp ! local value of psgn, double precision
INTEGER :: itmp ! temporary integer
CHARACTER(LEN=256) :: clinfo ! info character
CHARACTER(LEN=256) :: clname ! file name
CHARACTER(LEN=1) :: clrankpv, cldmspc !
CHARACTER(LEN=1) :: cl_type ! local value of cd_type
LOGICAL :: ll_only3rd ! T => if kstart, kcount present then *only* use values for 3rd spatial dimension.
+ LOGICAL :: llis1d, llis2d, llis3d
+ LOGICAL :: llsp ! use single precision
INTEGER :: inlev ! number of levels for 3D data
- REAL(dp) :: gma, gmi
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zwrk2d ! temporary arrays for reading into an array
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: zwrk3d ! smaller than jpi-jpj (XIOS only)
!---------------------------------------------------------------------
CHARACTER(LEN=lc) :: context
!
CALL set_xios_context(kiomid, context)
- ! Array shape information
- inlev = -1
- IF( PRESENT(pv_r1d) ) THEN
- irankpv = 1
- ishape(1:1) = SHAPE(pv_r1d)
- ELSE IF( PRESENT(pv_r2d) ) THEN
- irankpv = 2
- ishape(1:2) = SHAPE(pv_r2d)
- ELSE IF( PRESENT(pv_r3d) ) THEN
- irankpv = 3
- ishape(1:3) = SHAPE(pv_r3d)
- inlev = ishape(3)
+ !
+ llsp = PRESENT(pvsp1d) .OR. PRESENT(pvsp2d) .OR. PRESENT(pvsp3d)
+ IF( llsp ) THEN
+ llis1d = PRESENT(pvsp1d) ; IF( llis1d ) ishape(1:1) = SHAPE(pvsp1d)
+ llis2d = PRESENT(pvsp2d) ; IF( llis2d ) ishape(1:2) = SHAPE(pvsp2d)
+ llis3d = PRESENT(pvsp3d) ; IF( llis3d ) ishape(1:3) = SHAPE(pvsp3d)
+ ELSE
+ llis1d = PRESENT(pvdp1d) ; IF( llis1d ) ishape(1:1) = SHAPE(pvdp1d)
+ llis2d = PRESENT(pvdp2d) ; IF( llis2d ) ishape(1:2) = SHAPE(pvdp2d)
+ llis3d = PRESENT(pvdp3d) ; IF( llis3d ) ishape(1:3) = SHAPE(pvdp3d)
ENDIF
+ inlev = -1
+ IF( llis3d ) inlev = ishape(3)
!
idom = kdom
istop = nstop
@@ -1301,7 +1293,9 @@ CONTAINS
! else: if the file name finishes with _xxxx.nc with xxxx any number
ind1 = INDEX( clname, '_', back = .TRUE. ) + 1
ind2 = INDEX( clname, '.', back = .TRUE. ) - 1
- IF( ind2 > ind1 ) THEN ; IF( VERIFY( clname(ind1:ind2), '0123456789' ) == 0 ) idom = jpdom_local ; ENDIF
+ IF( ind2 > ind1 ) THEN
+ IF( VERIFY( clname(ind1:ind2), '0123456789' ) == 0 ) idom = jpdom_local
+ ENDIF
ENDIF
!
! check the consistency between input array and data rank in the file
@@ -1310,36 +1304,37 @@ CONTAINS
itime = 1
IF( PRESENT(ktime) ) itime = ktime
!
+ irankpv = 1 * COUNT( (/ llis1d /) ) + 2 * COUNT( (/ llis2d /) ) + 3 * COUNT( (/ llis3d /) )
WRITE(clrankpv, fmt='(i1)') irankpv
WRITE(cldmspc , fmt='(i1)') idmspc
!
IF( idmspc < irankpv ) THEN ! it seems we want to read more than we can...
- IF( irankpv == 3 .AND. idmspc == 2 ) THEN ! 3D input array from 2D spatial data in the file:
+ IF( llis3d .AND. idmspc == 2 ) THEN ! 3D input array from 2D spatial data in the file:
llok = inlev == 1 ! -> 3rd dimension must be equal to 1
- ELSEIF( irankpv == 3 .AND. idmspc == 1 ) THEN ! 3D input array from 1D spatial data in the file:
- llok = inlev == 1 .AND. SIZE(pv_r3d, 2) == 1 ! -> 2nd and 3rd dimensions must be equal to 1
- ELSEIF( irankpv == 2 .AND. idmspc == 2 ) THEN ! 2D input array from 1D spatial data in the file:
- llok = SIZE(pv_r2d, 2) == 1 ! -> 2nd dimension must be equal to 1
+ ELSEIF( llis3d .AND. idmspc == 1 ) THEN ! 3D input array from 1D spatial data in the file:
+ llok = inlev == 1 .AND. ishape(2) == 1 ! -> 2nd and 3rd dimensions must be equal to 1
+ ELSEIF( llis3d .AND. idmspc == 2 ) THEN ! 2D input array from 1D spatial data in the file:
+ llok = ishape(2) == 1 ! -> 2nd dimension must be equal to 1
ELSE
llok = .FALSE.
ENDIF
IF( .NOT. llok ) CALL ctl_stop( TRIM(clinfo), 'The file has only '//cldmspc//' spatial dimension', &
& '=> cannot read a true '//clrankpv//'D array from this file...' )
ELSEIF( idmspc == irankpv ) THEN
- IF( PRESENT(pv_r1d) .AND. idom /= jpdom_unknown ) &
+ IF( llis1d .AND. idom /= jpdom_unknown ) &
& CALL ctl_stop( TRIM(clinfo), 'case not coded...You must use jpdom_unknown' )
ELSEIF( idmspc > irankpv ) THEN ! it seems we want to read less than we should...
- IF( PRESENT(pv_r2d) .AND. itime == 1 .AND. idimsz(3) == 1 .AND. idmspc == 3 ) THEN
- CALL ctl_warn( TRIM(clinfo), '2D array input but 3 spatial dimensions in the file...' , &
- & 'As the size of the z dimension is 1 and as we try to read the first record, ', &
- & 'we accept this case, even if there is a possible mix-up between z and time dimension' )
- idmspc = idmspc - 1
+ IF( llis2d .AND. itime == 1 .AND. idimsz(3) == 1 .AND. idmspc == 3 ) THEN
+ CALL ctl_warn( TRIM(clinfo), '2D array input but 3 spatial dimensions in the file...' , &
+ & 'As the size of the z dimension is 1 and as we try to read the first record, ', &
+ & 'we accept this case, even if there is a possible mix-up between z and time dimension' )
+ idmspc = idmspc - 1
!!GS: possibility to read 3D ABL atmopsheric forcing and use 1st level to force BULK simulation
!ELSE
! CALL ctl_stop( TRIM(clinfo), 'To keep iom lisibility, when reading a '//clrankpv//'D array,', &
! & 'we do not accept data with '//cldmspc//' spatial dimensions' , &
! & 'Use ncwa -a to suppress the unnecessary dimensions' )
- ENDIF
+ ENDIF
ENDIF
!
! definition of istart and icnt
@@ -1355,11 +1350,11 @@ CONTAINS
ELSE
icnt (1:idmspc) = idimsz(1:idmspc)
ENDIF
- ELSE ! not a 1D array as pv_r1d requires jpdom_unknown
+ ELSE ! not a 1D array as pv(sd)p1d requires jpdom_unknown
! we do not read the overlap and the extra-halos -> from Nis0 to Nie0 and from Njs0 to Nje0
IF( idom == jpdom_global ) istart(1:2) = (/ mig(Nis0,0), mjg(Njs0,0) /)
icnt(1:2) = (/ Ni_0, Nj_0 /)
- IF( PRESENT(pv_r3d) ) THEN
+ IF( llis3d ) THEN
IF( idom == jpdom_auto_xy ) THEN
istart(3) = kstart(3)
icnt (3) = kcount(3)
@@ -1404,7 +1399,11 @@ CONTAINS
ENDIF
DO jl = 1, irankpv
WRITE( ctmp2, FMT="(', ', i1,'): ', i5,' /= icnt(', i1,'):', i5)" ) jl, ishape(jl), jl, icnt(jl)
- IF( ishape(jl) /= icnt(jl) ) CALL ctl_stop( TRIM(clinfo), 'size(pv_r'//clrankpv//TRIM(ctmp1)//TRIM(ctmp2) )
+ IF( llsp ) THEN
+ IF( ishape(jl) /= icnt(jl) ) CALL ctl_stop( TRIM(clinfo), 'size(pvsp'//clrankpv//TRIM(ctmp1)//TRIM(ctmp2) )
+ ELSE
+ IF( ishape(jl) /= icnt(jl) ) CALL ctl_stop( TRIM(clinfo), 'size(pvdp'//clrankpv//TRIM(ctmp1)//TRIM(ctmp2) )
+ ENDIF
END DO
ENDIF
@@ -1413,21 +1412,30 @@ CONTAINS
!-
IF( idvar > 0 .AND. istop == nstop ) THEN ! no additional errors until this point...
!
- CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, pv_r3d )
+ CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, &
+ & pvsp1d, pvsp2d, pvsp3d, pvdp1d, pvdp2d, pvdp3d )
IF( istop == nstop ) THEN ! no additional errors until this point...
IF(lwp) WRITE(numout,"(10x,' read ',a,' (rec: ',i6,') in ',a,' ok')") TRIM(cdvar), itime, TRIM(iom_file(kiomid)%name)
+ !--- overlap areas and extra hallows (mpp)
cl_type = 'T'
IF( PRESENT(cd_type) ) cl_type = cd_type
- !--- halos and NP folding (NP folding to be done even if no halos)
- IF( idom /= jpdom_unknown .AND. cl_type /= 'Z' .AND. ( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) ) THEN
- zsgn = 1._wp
- IF( PRESENT(psgn ) ) zsgn = psgn
- IF( PRESENT(pv_r2d) ) THEN
- CALL lbc_lnk( 'iom', pv_r2d, cl_type, zsgn, kfillmode = kfill, ldfull = .TRUE. )
- ELSEIF( PRESENT(pv_r3d) ) THEN
- CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill, ldfull = .TRUE. )
+ IF( llsp ) THEN
+ zsgn_sp = 1._sp
+ IF( PRESENT(psgn_sp) ) zsgn_sp = psgn_sp
+ IF( llis2d .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+ CALL lbc_lnk( 'iom', pvsp2d, cl_type, zsgn_sp, kfillmode = kfill, ldfull = .TRUE. )
+ ELSEIF( llis3d .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+ CALL lbc_lnk( 'iom', pvsp3d, cl_type, zsgn_sp, kfillmode = kfill, ldfull = .TRUE. )
+ ENDIF
+ ELSE
+ zsgn_dp = 1._dp
+ IF( PRESENT(psgn_dp) ) zsgn_dp = psgn_dp
+ IF( llis2d .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+ CALL lbc_lnk( 'iom', pvdp2d, cl_type, zsgn_dp, kfillmode = kfill, ldfull = .TRUE. )
+ ELSEIF( llis3d .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+ CALL lbc_lnk( 'iom', pvdp3d, cl_type, zsgn_dp, kfillmode = kfill, ldfull = .TRUE. )
ENDIF
ENDIF
!
@@ -1442,52 +1450,69 @@ CONTAINS
!
ELSE ! read using XIOS. Only if key_xios is defined
#if defined key_xios
-!would be good to be able to check which context is active and swap only if current is not restart
+ !would be good to be able to check which context is active and swap only if current is not restart
idvar = iom_varid( kiomid, cdvar )
CALL iom_swap(context)
- zsgn = 1._wp
- IF( PRESENT(psgn ) ) zsgn = psgn
+ IF( llsp ) THEN
+ zsgn_sp = 1._sp ; IF( PRESENT(psgn_sp) ) zsgn_sp = psgn_sp
+ ELSE
+ zsgn_dp = 1._dp ; IF( PRESENT(psgn_dp) ) zsgn_dp = psgn_dp
+ ENDIF
cl_type = 'T'
IF( PRESENT(cd_type) ) cl_type = cd_type
- ! We do not know whether pv_r[23]d will have halo points during context initialisation in iom_set_vars_active,
- ! so we must always read to an array on the full grid (1:jpi,1:jpj), then copy the correct part to pv_r[23]d
+ ! We do not know whether pv[sd]p[23]d will have halo points during context initialisation in iom_set_vars_active,
+ ! so we must always read to an array on the full grid (1:jpi,1:jpj), then copy the correct part to pv[sd]p[23]d
IF( irankpv > 1 .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
- llok = (ishape(1) == jpi) .AND. (ishape(2) == jpj)
+ llok = ishape(1) == jpi .AND. ishape(2) == jpj
ihls = arr_hls( ishape(1), ishape(2), ldtile=.FALSE., ldsize=.TRUE. )
ELSE
llok = .TRUE.
ENDIF
- IF( PRESENT(pv_r3d) ) THEN
+ IF( llis3d ) THEN
IF(lwp) WRITE(numout,*) 'XIOS RST READ (3D): ',TRIM(cdvar)
IF( .NOT. llok ) THEN
ALLOCATE( zwrk3d(jpi,jpj,inlev) )
CALL xios_recv_field( trim(cdvar), zwrk3d(:,:,:) )
- pv_r3d(:,:,:) = zwrk3d(A1Di(ihls(1)),A1Dj(ihls(2)),:)
+ IF( llsp ) THEN ; pvsp3d(:,:,:) = zwrk3d(A1Di(ihls(1)),A1Dj(ihls(2)),:)
+ ELSE ; pvdp3d(:,:,:) = zwrk3d(A1Di(ihls(1)),A1Dj(ihls(2)),:)
+ ENDIF
DEALLOCATE( zwrk3d )
ELSE
- CALL xios_recv_field( trim(cdvar), pv_r3d(:,:,:) )
+ IF( llsp ) THEN ; CALL xios_recv_field( TRIM(cdvar), pvsp3d )
+ ELSE ; CALL xios_recv_field( TRIM(cdvar), pvdp3d )
+ ENDIF
ENDIF
IF(idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
- CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill, ldfull = .TRUE. )
+ IF( llsp ) THEN ; CALL lbc_lnk( 'iom', pvsp3d, cl_type, zsgn_sp, kfillmode = kfill, ldfull = .TRUE. )
+ ELSE ; CALL lbc_lnk( 'iom', pvdp3d, cl_type, zsgn_dp, kfillmode = kfill, ldfull = .TRUE. )
+ ENDIF
ENDIF
- ELSEIF( PRESENT(pv_r2d) ) THEN
+ ELSEIF( llis2d ) THEN
IF(lwp) WRITE(numout,*) 'XIOS RST READ (2D): ', TRIM(cdvar)
IF( .NOT. llok ) THEN
ALLOCATE( zwrk2d(jpi,jpj) )
- CALL xios_recv_field( trim(cdvar), zwrk2d(:,:) )
- pv_r2d(:,:) = zwrk2d(A1Di(ihls(1)),A1Dj(ihls(2)))
+ CALL xios_recv_field( TRIM(cdvar), zwrk2d(:,:) )
+ IF( llsp ) THEN ; pvsp2d(:,:) = zwrk2d(A1Di(ihls(1)),A1Dj(ihls(2)))
+ ELSE ; pvdp2d(:,:) = zwrk2d(A1Di(ihls(1)),A1Dj(ihls(2)))
+ ENDIF
DEALLOCATE( zwrk2d )
ELSE
- CALL xios_recv_field( trim(cdvar), pv_r2d(:,:) )
+ IF( llsp ) THEN ; CALL xios_recv_field( TRIM(cdvar), pvsp2d )
+ ELSE ; CALL xios_recv_field( TRIM(cdvar), pvdp2d )
+ ENDIF
ENDIF
IF(idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
- CALL lbc_lnk('iom', pv_r2d, cl_type, zsgn, kfillmode = kfill, ldfull = .TRUE. )
+ IF( llsp ) THEN ; CALL lbc_lnk('iom', pvsp2d, cl_type, zsgn_sp, kfillmode = kfill, ldfull = .TRUE. )
+ ELSE ; CALL lbc_lnk('iom', pvdp2d, cl_type, zsgn_dp, kfillmode = kfill, ldfull = .TRUE. )
+ ENDIF
ENDIF
- ELSEIF( PRESENT(pv_r1d) ) THEN
+ ELSEIF( llis1d ) THEN
IF(lwp) WRITE(numout,*) 'XIOS RST READ (1D): ', TRIM(cdvar)
- CALL xios_recv_field( TRIM(cdvar), pv_r1d)
+ IF( llsp ) THEN ; CALL xios_recv_field( TRIM(cdvar), pvsp1d )
+ ELSE ; CALL xios_recv_field( TRIM(cdvar), pvdp1d )
+ ENDIF
ENDIF
CALL iom_swap(cxios_context)
#else
@@ -1497,17 +1522,32 @@ CONTAINS
ENDIF
!--- Apply scale_factor and offset
- zscf = iom_file(kiomid)%scf(idvar) ! scale factor
- zofs = iom_file(kiomid)%ofs(idvar) ! offset
- IF( PRESENT(pv_r1d) ) THEN
- IF( zscf /= 1._wp ) pv_r1d(:) = pv_r1d(:) * zscf
- IF( zofs /= 0._wp ) pv_r1d(:) = pv_r1d(:) + zofs
- ELSEIF( PRESENT(pv_r2d) ) THEN
- IF( zscf /= 1._wp) pv_r2d(:,:) = pv_r2d(:,:) * zscf
- IF( zofs /= 0._wp) pv_r2d(:,:) = pv_r2d(:,:) + zofs
- ELSEIF( PRESENT(pv_r3d) ) THEN
- IF( zscf /= 1._wp) pv_r3d(:,:,:) = pv_r3d(:,:,:) * zscf
- IF( zofs /= 0._wp) pv_r3d(:,:,:) = pv_r3d(:,:,:) + zofs
+ IF( llsp ) THEN
+ zscf_sp = iom_file(kiomid)%scf(idvar) ! scale factor
+ zofs_sp = iom_file(kiomid)%ofs(idvar) ! offset
+ IF( llis1d ) THEN
+ IF( zscf_sp /= 1._sp ) pvsp1d(: ) = pvsp1d(: ) * zscf_sp
+ IF( zofs_sp /= 0._sp ) pvsp1d(: ) = pvsp1d(: ) + zofs_sp
+ ELSEIF( llis2d ) THEN
+ IF( zscf_sp /= 1._sp ) pvsp2d(:,: ) = pvsp2d(:,: ) * zscf_sp
+ IF( zofs_sp /= 0._sp ) pvsp2d(:,: ) = pvsp2d(:,: ) + zofs_sp
+ ELSEIF( llis3d ) THEN
+ IF( zscf_sp /= 1._sp ) pvsp3d(:,:,:) = pvsp3d(:,:,:) * zscf_sp
+ IF( zofs_sp /= 0._sp ) pvsp3d(:,:,:) = pvsp3d(:,:,:) + zofs_sp
+ ENDIF
+ ELSE
+ zscf_dp = iom_file(kiomid)%scf(idvar) ! scale factor
+ zofs_dp = iom_file(kiomid)%ofs(idvar) ! offset
+ IF( llis1d ) THEN
+ IF( zscf_dp /= 1._dp ) pvdp1d(: ) = pvdp1d(: ) * zscf_dp
+ IF( zofs_dp /= 0._dp ) pvdp1d(: ) = pvdp1d(: ) + zofs_dp
+ ELSEIF( llis2d ) THEN
+ IF( zscf_dp /= 1._dp ) pvdp2d(:,: ) = pvdp2d(:,: ) * zscf_dp
+ IF( zofs_dp /= 0._dp ) pvdp2d(:,: ) = pvdp2d(:,: ) + zofs_dp
+ ELSEIF( llis3d ) THEN
+ IF( zscf_dp /= 1._dp ) pvdp3d(:,:,:) = pvdp3d(:,:,:) * zscf_dp
+ IF( zofs_dp /= 0._dp ) pvdp3d(:,:,:) = pvdp3d(:,:,:) + zofs_dp
+ ENDIF
ENDIF
!
END SUBROUTINE iom_get_123d
@@ -1693,32 +1733,8 @@ CONTAINS
INTEGER :: ivid ! variable id
CHARACTER(LEN=lc) :: context
!
- CALL set_xios_context(kiomid, context)
-
- llx = .NOT. (context == "NONE")
-
- IF( llx ) THEN
-#ifdef key_xios
- IF( kt == kwrite ) THEN
- IF(lwp) write(numout,*) 'RESTART: write (XIOS 0D) ',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_put(TRIM(cdvar), pvar)
- CALL iom_swap(cxios_context)
- ELSE
- IF(lwp) write(numout,*) 'RESTART: define (XIOS 0D) ',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_set_rstw_active( TRIM(cdvar), rs0 = pvar )
- CALL iom_swap(cxios_context)
- ENDIF
-#endif
- ELSE
- IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) THEN
- ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
- CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r0d = real(pvar, dp) )
- ENDIF
- ENDIF
- ENDIF
+ CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvsp0d = pvar )
+ !
END SUBROUTINE iom_rp0d_sp
SUBROUTINE iom_rp0d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
@@ -1733,32 +1749,8 @@ CONTAINS
INTEGER :: ivid ! variable id
CHARACTER(LEN=lc) :: context
!
- CALL set_xios_context(kiomid, context)
-
- llx = .NOT. (context == "NONE")
-
- IF( llx ) THEN
-#ifdef key_xios
- IF( kt == kwrite ) THEN
- IF(lwp) write(numout,*) 'RESTART: write (XIOS 0D) ',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_put(TRIM(cdvar), pvar)
- CALL iom_swap(cxios_context)
- ELSE
- IF(lwp) write(numout,*) 'RESTART: define (XIOS 0D) ',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_set_rstw_active( TRIM(cdvar), rd0 = pvar )
- CALL iom_swap(cxios_context)
- ENDIF
-#endif
- ELSE
- IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) THEN
- ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
- CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r0d = pvar )
- ENDIF
- ENDIF
- ENDIF
+ CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvdp0d = pvar )
+ !
END SUBROUTINE iom_rp0d_dp
@@ -1774,32 +1766,8 @@ CONTAINS
INTEGER :: ivid ! variable id
CHARACTER(LEN=lc) :: context
!
- CALL set_xios_context(kiomid, context)
-
- llx = .NOT. (context == "NONE")
-
- IF( llx ) THEN
-#ifdef key_xios
- IF( kt == kwrite ) THEN
- IF(lwp) write(numout,*) 'RESTART: write (XIOS 1D) ',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_put(TRIM(cdvar), pvar)
- CALL iom_swap(cxios_context)
- ELSE
- IF(lwp) write(numout,*) 'RESTART: define (XIOS 1D)',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_set_rstw_active( TRIM(cdvar), rs1 = pvar )
- CALL iom_swap(cxios_context)
- ENDIF
-#endif
- ELSE
- IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) THEN
- ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
- CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r1d = real(pvar, dp) )
- ENDIF
- ENDIF
- ENDIF
+ CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvsp1d = pvar )
+ !
END SUBROUTINE iom_rp1d_sp
SUBROUTINE iom_rp1d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
@@ -1814,35 +1782,11 @@ CONTAINS
INTEGER :: ivid ! variable id
CHARACTER(LEN=lc) :: context
!
- CALL set_xios_context(kiomid, context)
-
- llx = .NOT. (context == "NONE")
-
- IF( llx ) THEN
-#ifdef key_xios
- IF( kt == kwrite ) THEN
- IF(lwp) write(numout,*) 'RESTART: write (XIOS 1D) ',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_put(TRIM(cdvar), pvar)
- CALL iom_swap(cxios_context)
- ELSE
- IF(lwp) write(numout,*) 'RESTART: define (XIOS 1D)',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_set_rstw_active( TRIM(cdvar), rd1 = pvar )
- CALL iom_swap(cxios_context)
- ENDIF
-#endif
- ELSE
- IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) THEN
- ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
- CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r1d = pvar )
- ENDIF
- ENDIF
- ENDIF
+ CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvdp1d = pvar )
+ !
END SUBROUTINE iom_rp1d_dp
-
+
SUBROUTINE iom_rp2d_sp( kt, kwrite, kiomid, cdvar, pvar, ktype )
INTEGER , INTENT(in) :: kt ! ocean time-step
INTEGER , INTENT(in) :: kwrite ! writing time-step
@@ -1855,32 +1799,8 @@ CONTAINS
INTEGER :: ivid ! variable id
CHARACTER(LEN=lc) :: context
!
- CALL set_xios_context(kiomid, context)
-
- llx = .NOT. (context == "NONE")
-
- IF( llx ) THEN
-#ifdef key_xios
- IF( kt == kwrite ) THEN
- IF(lwp) write(numout,*) 'RESTART: write (XIOS 2D) ',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_put(TRIM(cdvar), pvar)
- CALL iom_swap(cxios_context)
- ELSE
- IF(lwp) write(numout,*) 'RESTART: define (XIOS 2D)',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_set_rstw_active( TRIM(cdvar), rs2 = pvar )
- CALL iom_swap(cxios_context)
- ENDIF
-#endif
- ELSE
- IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) THEN
- ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
- CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r2d = real(pvar, dp) )
- ENDIF
- ENDIF
- ENDIF
+ CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvsp2d = pvar )
+ !
END SUBROUTINE iom_rp2d_sp
SUBROUTINE iom_rp2d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
@@ -1895,32 +1815,8 @@ CONTAINS
INTEGER :: ivid ! variable id
CHARACTER(LEN=lc) :: context
!
- CALL set_xios_context(kiomid, context)
-
- llx = .NOT. (context == "NONE")
-
- IF( llx ) THEN
-#ifdef key_xios
- IF( kt == kwrite ) THEN
- IF(lwp) write(numout,*) 'RESTART: write (XIOS 2D) ',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_put(TRIM(cdvar), pvar)
- CALL iom_swap(cxios_context)
- ELSE
- IF(lwp) write(numout,*) 'RESTART: define (XIOS 2D)',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_set_rstw_active( TRIM(cdvar), rd2 = pvar )
- CALL iom_swap(cxios_context)
- ENDIF
-#endif
- ELSE
- IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) THEN
- ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
- CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r2d = pvar )
- ENDIF
- ENDIF
- ENDIF
+ CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvdp2d = pvar )
+ !
END SUBROUTINE iom_rp2d_dp
@@ -1936,32 +1832,8 @@ CONTAINS
INTEGER :: ivid ! variable id
CHARACTER(LEN=lc) :: context
!
- CALL set_xios_context(kiomid, context)
-
- llx = .NOT. (context == "NONE")
-
- IF( llx ) THEN
-#ifdef key_xios
- IF( kt == kwrite ) THEN
- IF(lwp) write(numout,*) 'RESTART: write (XIOS 3D) ',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_put(TRIM(cdvar), pvar)
- CALL iom_swap(cxios_context)
- ELSE
- IF(lwp) write(numout,*) 'RESTART: define (XIOS 3D)',TRIM(cdvar)
- CALL iom_swap(context)
- CALL iom_set_rstw_active( TRIM(cdvar), rs3 = pvar )
- CALL iom_swap(cxios_context)
- ENDIF
-#endif
- ELSE
- IF( kiomid > 0 ) THEN
- IF( iom_file(kiomid)%nfid > 0 ) THEN
- ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
- CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r3d = real(pvar, dp) )
- ENDIF
- ENDIF
- ENDIF
+ CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvsp3d = pvar )
+ !
END SUBROUTINE iom_rp3d_sp
SUBROUTINE iom_rp3d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
@@ -1976,6 +1848,29 @@ CONTAINS
INTEGER :: ivid ! variable id
CHARACTER(LEN=lc) :: context
!
+ CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvdp3d = pvar )
+ !
+ END SUBROUTINE iom_rp3d_dp
+
+ SUBROUTINE iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvsp0d, pvsp1d, pvsp2d, pvsp3d, pvdp0d, pvdp1d, pvdp2d, pvdp3d )
+ INTEGER , INTENT(in) :: kt ! ocean time-step
+ INTEGER , INTENT(in) :: kwrite ! writing time-step
+ INTEGER , INTENT(in) :: kiomid ! Identifier of the file
+ CHARACTER(len=*) , INTENT(in) :: cdvar ! time axis name
+ INTEGER , INTENT(in), OPTIONAL :: ktype ! variable external type
+ REAL(sp) , INTENT(in), OPTIONAL :: pvsp0d ! read field (0D case), single precision
+ REAL(sp), DIMENSION(:) , INTENT(in), OPTIONAL :: pvsp1d ! read field (1D case), single precision
+ REAL(sp), DIMENSION(:,:) , INTENT(in), OPTIONAL :: pvsp2d ! read field (2D case), single precision
+ REAL(sp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: pvsp3d ! read field (3D case), single precision
+ REAL(dp) , INTENT(in), OPTIONAL :: pvdp0d ! read field (0D case), double precision
+ REAL(dp), DIMENSION(:) , INTENT(in), OPTIONAL :: pvdp1d ! read field (1D case), double precision
+ REAL(dp), DIMENSION(:,:) , INTENT(in), OPTIONAL :: pvdp2d ! read field (2D case), double precision
+ REAL(dp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL :: pvdp3d ! read field (3D case), double precision
+ !
+ LOGICAL :: llx ! local xios write flag
+ INTEGER :: ivid ! variable id
+ CHARACTER(LEN=lc) :: context
+ !
CALL set_xios_context(kiomid, context)
llx = .NOT. (context == "NONE")
@@ -1983,14 +1878,22 @@ CONTAINS
IF( llx ) THEN
#ifdef key_xios
IF( kt == kwrite ) THEN
- IF(lwp) write(numout,*) 'RESTART: write (XIOS 3D) ',TRIM(cdvar)
+ IF(lwp) write(numout,*) 'RESTART: write (XIOS) ',TRIM(cdvar)
CALL iom_swap(context)
- CALL iom_put(TRIM(cdvar), pvar)
+ IF( PRESENT(pvsp0d) ) CALL iom_put(TRIM(cdvar), pvsp0d)
+ IF( PRESENT(pvsp1d) ) CALL iom_put(TRIM(cdvar), pvsp1d)
+ IF( PRESENT(pvsp2d) ) CALL iom_put(TRIM(cdvar), pvsp2d)
+ IF( PRESENT(pvsp3d) ) CALL iom_put(TRIM(cdvar), pvsp3d)
+ IF( PRESENT(pvdp0d) ) CALL iom_put(TRIM(cdvar), pvdp0d)
+ IF( PRESENT(pvdp1d) ) CALL iom_put(TRIM(cdvar), pvdp1d)
+ IF( PRESENT(pvdp2d) ) CALL iom_put(TRIM(cdvar), pvdp2d)
+ IF( PRESENT(pvdp3d) ) CALL iom_put(TRIM(cdvar), pvdp3d)
CALL iom_swap(cxios_context)
ELSE
- IF(lwp) write(numout,*) 'RESTART: define (XIOS 3D)',TRIM(cdvar)
+ IF(lwp) write(numout,*) 'RESTART: define (XIOS)',TRIM(cdvar)
CALL iom_swap(context)
- CALL iom_set_rstw_active( TRIM(cdvar), rd3 = pvar )
+ CALL iom_set_rstw_active( TRIM(cdvar), rs0 = pvsp0d, rs1 = pvsp1d, rs2 = pvsp2d, rs3 = pvsp3d &
+ & , rd0 = pvdp0d, rd1 = pvdp1d, rd2 = pvdp2d, rd3 = pvdp3d )
CALL iom_swap(cxios_context)
ENDIF
#endif
@@ -1998,12 +1901,12 @@ CONTAINS
IF( kiomid > 0 ) THEN
IF( iom_file(kiomid)%nfid > 0 ) THEN
ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
- CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r3d = pvar )
+ CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pvsp0d, pvsp1d, pvsp2d, pvsp3d, &
+ & pvdp0d, pvdp1d, pvdp2d, pvdp3d )
ENDIF
ENDIF
ENDIF
- END SUBROUTINE iom_rp3d_dp
-
+ END SUBROUTINE iom_rp0123d
SUBROUTINE iom_delay_rst( cdaction, cdcpnt, kncid )
diff --git a/src/OCE/IOM/iom_nf90.F90 b/src/OCE/IOM/iom_nf90.F90
index 38d9943cd6de362e2e7beb650e2b0de87b7e2e57..a88ad855a517e68ea016dbd8b5962cf01fa58c0d 100644
--- a/src/OCE/IOM/iom_nf90.F90
+++ b/src/OCE/IOM/iom_nf90.F90
@@ -34,10 +34,7 @@ MODULE iom_nf90
INTERFACE iom_nf90_get
MODULE PROCEDURE iom_nf90_g0d_sp
- MODULE PROCEDURE iom_nf90_g0d_dp, iom_nf90_g123d_dp
- END INTERFACE
- INTERFACE iom_nf90_rstput
- MODULE PROCEDURE iom_nf90_rp0123d_dp
+ MODULE PROCEDURE iom_nf90_g0d_dp, iom_nf90_g123d
END INTERFACE
!! * Substitutions
@@ -312,8 +309,8 @@ CONTAINS
CALL iom_nf90_check(NF90_GET_VAR(iom_file(kiomid)%nfid, iom_file(kiomid)%nvid(kvid), pvar, start = kstart), clinfo )
END SUBROUTINE iom_nf90_g0d_dp
- SUBROUTINE iom_nf90_g123d_dp( kiomid, kvid, knbdim, kstart, kcount, kx1, kx2, ky1, ky2, &
- & pv_r1d, pv_r2d, pv_r3d )
+ SUBROUTINE iom_nf90_g123d( kiomid, kvid, knbdim, kstart, kcount, kx1, kx2, ky1, ky2, &
+ & pvsp1d, pvsp2d, pvsp3d, pvdp1d, pvdp2d, pvdp3d )
!!-----------------------------------------------------------------------
!! *** ROUTINE iom_nf90_g123d ***
!!
@@ -327,9 +324,12 @@ CONTAINS
INTEGER , DIMENSION(:) , INTENT(in ) :: kstart ! start position of the reading in each axis
INTEGER , DIMENSION(:) , INTENT(in ) :: kcount ! number of points to be read in each axis
INTEGER , INTENT(in ) :: kx1, kx2, ky1, ky2 ! subdomain indexes
- REAL(dp), DIMENSION(:) , INTENT( out), OPTIONAL :: pv_r1d ! read field (1D case)
- REAL(dp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pv_r2d ! read field (2D case)
- REAL(dp), DIMENSION(:,:,:) , INTENT( out), OPTIONAL :: pv_r3d ! read field (3D case)
+ REAL(sp), DIMENSION(:) , INTENT( out), OPTIONAL :: pvsp1d ! read field (1D case), single precision
+ REAL(sp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pvsp2d ! read field (2D case), single precision
+ REAL(sp), DIMENSION(:,:,:) , INTENT( out), OPTIONAL :: pvsp3d ! read field (3D case), single precision
+ REAL(dp), DIMENSION(:) , INTENT( out), OPTIONAL :: pvdp1d ! read field (1D case), double precision
+ REAL(dp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pvdp2d ! read field (2D case), double precision
+ REAL(dp), DIMENSION(:,:,:) , INTENT( out), OPTIONAL :: pvdp3d ! read field (3D case), double precision
!
CHARACTER(LEN=100) :: clinfo ! info character
INTEGER :: if90id ! nf90 identifier of the opened file
@@ -339,18 +339,27 @@ CONTAINS
if90id = iom_file(kiomid)%nfid ! get back NetCDF file id
ivid = iom_file(kiomid)%nvid(kvid) ! get back NetCDF var id
!
- IF( PRESENT(pv_r1d) ) THEN
- CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pv_r1d(: ), start = kstart(1:knbdim), &
+ IF( PRESENT(pvsp1d) ) THEN
+ CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvsp1d(: ), start = kstart(1:knbdim), &
+ & count = kcount(1:knbdim)), clinfo )
+ ELSEIF( PRESENT(pvdp1d) ) THEN
+ CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvdp1d(: ), start = kstart(1:knbdim), &
+ & count = kcount(1:knbdim)), clinfo )
+ ELSEIF( PRESENT(pvsp2d) ) THEN
+ CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvsp2d(kx1:kx2,ky1:ky2 ), start = kstart(1:knbdim), &
+ & count = kcount(1:knbdim)), clinfo )
+ ELSEIF( PRESENT(pvdp2d) ) THEN
+ CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvdp2d(kx1:kx2,ky1:ky2 ), start = kstart(1:knbdim), &
& count = kcount(1:knbdim)), clinfo )
- ELSEIF( PRESENT(pv_r2d) ) THEN
- CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pv_r2d(kx1:kx2,ky1:ky2 ), start = kstart(1:knbdim), &
+ ELSEIF( PRESENT(pvsp3d) ) THEN
+ CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvsp3d(kx1:kx2,ky1:ky2,:), start = kstart(1:knbdim), &
& count = kcount(1:knbdim)), clinfo )
- ELSEIF( PRESENT(pv_r3d) ) THEN
- CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pv_r3d(kx1:kx2,ky1:ky2,:), start = kstart(1:knbdim), &
+ ELSEIF( PRESENT(pvdp3d) ) THEN
+ CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvdp3d(kx1:kx2,ky1:ky2,:), start = kstart(1:knbdim), &
& count = kcount(1:knbdim)), clinfo )
ENDIF
!
- END SUBROUTINE iom_nf90_g123d_dp
+ END SUBROUTINE iom_nf90_g123d
@@ -524,8 +533,9 @@ CONTAINS
!
END SUBROUTINE iom_nf90_putatt
- SUBROUTINE iom_nf90_rp0123d_dp( kt, kwrite, kiomid, cdvar , kvid , ktype, &
- & pv_r0d, pv_r1d, pv_r2d, pv_r3d )
+ SUBROUTINE iom_nf90_rstput( kt, kwrite, kiomid, cdvar , kvid , ktype , &
+ & pvsp0d, pvsp1d, pvsp2d, pvsp3d, &
+ & pvdp0d, pvdp1d, pvdp2d, pvdp3d )
!!--------------------------------------------------------------------
!! *** SUBROUTINE iom_nf90_rstput ***
!!
@@ -537,10 +547,14 @@ CONTAINS
CHARACTER(len=*) , INTENT(in) :: cdvar ! variable name
INTEGER , INTENT(in) :: kvid ! variable id
INTEGER , INTENT(in), OPTIONAL :: ktype ! variable type (default R8)
- REAL(dp) , INTENT(in), OPTIONAL :: pv_r0d ! written Od field
- REAL(dp), DIMENSION( :), INTENT(in), OPTIONAL :: pv_r1d ! written 1d field
- REAL(dp), DIMENSION(:, : ), INTENT(in), OPTIONAL :: pv_r2d ! written 2d field
- REAL(dp), DIMENSION(:, :, :), INTENT(in), OPTIONAL :: pv_r3d ! written 3d field
+ REAL(sp) , INTENT(in), OPTIONAL :: pvsp0d ! written Od field
+ REAL(sp), DIMENSION( :), INTENT(in), OPTIONAL :: pvsp1d ! written 1d field
+ REAL(sp), DIMENSION(:, : ), INTENT(in), OPTIONAL :: pvsp2d ! written 2d field
+ REAL(sp), DIMENSION(:, :, :), INTENT(in), OPTIONAL :: pvsp3d ! written 3d field
+ REAL(dp) , INTENT(in), OPTIONAL :: pvdp0d ! written Od field
+ REAL(dp), DIMENSION( :), INTENT(in), OPTIONAL :: pvdp1d ! written 1d field
+ REAL(dp), DIMENSION(:, : ), INTENT(in), OPTIONAL :: pvdp2d ! written 2d field
+ REAL(dp), DIMENSION(:, :, :), INTENT(in), OPTIONAL :: pvdp3d ! written 3d field
!
INTEGER :: idims ! number of dimension
INTEGER :: idvar ! variable id
@@ -555,12 +569,18 @@ CONTAINS
! ! nn_nchunks_[i,j,k,t] namelist parameters
INTEGER :: ichunkalg, ishuffle, ideflate, ideflate_level
! ! NetCDF4 internally fixed parameters
+ INTEGER :: idlv ! local variable
LOGICAL :: lchunk ! logical switch to activate chunking and compression
! ! when appropriate (currently chunking is applied to 4d fields only)
- INTEGER :: idlv ! local variable
+ LOGICAL :: llis0d, llis1d, llis2d, llis3d
CHARACTER(LEN=256) :: ccname ! local variable
!---------------------------------------------------------------------
!
+ llis0d = PRESENT(pvsp0d) .OR. PRESENT(pvdp0d)
+ llis1d = PRESENT(pvsp1d) .OR. PRESENT(pvdp1d)
+ llis2d = PRESENT(pvsp2d) .OR. PRESENT(pvdp2d)
+ llis3d = PRESENT(pvsp3d) .OR. PRESENT(pvdp3d)
+ !
clinfo = ' iom_nf90_rp0123d, file: '//TRIM(iom_file(kiomid)%name)//', var: '//TRIM(cdvar)
if90id = iom_file(kiomid)%nfid
!
@@ -609,12 +629,12 @@ CONTAINS
CALL iom_nf90_check(NF90_REDEF( if90id ), clinfo) ; iom_file(kiomid)%irec = -1
ENDIF
! variable definition
- IF( PRESENT(pv_r0d) ) THEN ; idims = 0
- ELSEIF( PRESENT(pv_r1d) ) THEN
- idims = 2 ; idimid(1:idims) = (/3,4/)
- ELSEIF( PRESENT(pv_r2d) ) THEN ; idims = 3 ; idimid(1:idims) = (/1,2,4/)
- ELSEIF( PRESENT(pv_r3d) ) THEN
- idims = 4 ; idimid(1:idims) = (/1,2,3,4/)
+ IF( llis0d ) THEN ; idims = 0
+ ELSEIF( llis1d ) THEN
+ idims = 2 ; idimid(1:idims) = (/3,4/)
+ ELSEIF( llis2d ) THEN ; idims = 3 ; idimid(1:idims) = (/1,2,4/)
+ ELSEIF( llis3d ) THEN
+ idims = 4 ; idimid(1:idims) = (/1,2,3,4/)
ENDIF
IF( PRESENT(ktype) ) THEN ! variable external type
SELECT CASE (ktype)
@@ -628,7 +648,7 @@ CONTAINS
ELSE
itype = NF90_DOUBLE
ENDIF
- IF( PRESENT(pv_r0d) ) THEN
+ IF( llis0d ) THEN
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, iom_file(kiomid)%nvid(idvar) ), clinfo )
ELSE
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims), iom_file(kiomid)%nvid(idvar) ), clinfo )
@@ -641,9 +661,7 @@ CONTAINS
iom_file(kiomid)%scf(idvar) = 1.
iom_file(kiomid)%ofs(idvar) = 0.
iom_file(kiomid)%ndims(idvar) = idims
- IF( .NOT. PRESENT(pv_r0d) ) THEN ; iom_file(kiomid)%luld(idvar) = .TRUE.
- ELSE ; iom_file(kiomid)%luld(idvar) = .FALSE.
- ENDIF
+ iom_file(kiomid)%luld(idvar) = .NOT. llis0d
DO jd = 1, idims
CALL iom_nf90_check(NF90_INQUIRE_DIMENSION( if90id, idimid(jd), len = iom_file(kiomid)%dimsz(jd,idvar) ), clinfo)
IF ( lchunk ) ichunksz(jd) = iom_file(kiomid)%dimsz(jd,idvar)
@@ -672,7 +690,7 @@ CONTAINS
CALL iom_nf90_check(NF90_ENDDEF( if90id ), clinfo) ; iom_file(kiomid)%irec = 0
ENDIF
! on what kind of domain must the data be written?
- IF( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) THEN
+ IF( llis2d .OR. llis3d ) THEN
! write dimension variables if it is not already done
! =============
@@ -695,8 +713,10 @@ CONTAINS
IF(lwp) WRITE(numout,*) TRIM(clinfo)//' write dimension variables done'
ENDIF
- IF( PRESENT(pv_r2d) ) ishape(1:2) = SHAPE(pv_r2d)
- IF( PRESENT(pv_r3d) ) ishape(1:3) = SHAPE(pv_r3d)
+ IF( PRESENT(pvsp2d) ) ishape(1:2) = SHAPE(pvsp2d)
+ IF( PRESENT(pvdp2d) ) ishape(1:2) = SHAPE(pvdp2d)
+ IF( PRESENT(pvsp3d) ) ishape(1:3) = SHAPE(pvsp3d)
+ IF( PRESENT(pvdp3d) ) ishape(1:3) = SHAPE(pvdp3d)
IF( ishape(1) == Ni_0 .AND. ishape(2) == Nj_0 ) THEN ! array with 0 halo
ix1 = 1 ; ix2 = Ni_0 ; iy1 = 1 ; iy2 = Nj_0
ELSEIF( ishape(1) == jpi .AND. ishape(2) == jpj ) THEN ! array with nn_hls halos
@@ -711,14 +731,22 @@ CONTAINS
! write the data
! =============
- IF( PRESENT(pv_r0d) ) THEN
- CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r0d ), clinfo )
- ELSEIF( PRESENT(pv_r1d) ) THEN
- CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r1d(:) ), clinfo )
- ELSEIF( PRESENT(pv_r2d) ) THEN
- CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r2d(ix1:ix2,iy1:iy2) ), clinfo )
- ELSEIF( PRESENT(pv_r3d) ) THEN
- CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r3d(ix1:ix2,iy1:iy2,:) ), clinfo )
+ IF( PRESENT(pvsp0d) ) THEN
+ CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvsp0d ), clinfo )
+ ELSEIF( PRESENT(pvdp0d) ) THEN
+ CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvdp0d ), clinfo )
+ ELSEIF( PRESENT(pvsp1d) ) THEN
+ CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvsp1d(:) ), clinfo )
+ ELSEIF( PRESENT(pvdp1d) ) THEN
+ CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvdp1d(:) ), clinfo )
+ ELSEIF( PRESENT(pvsp2d) ) THEN
+ CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvsp2d(ix1:ix2,iy1:iy2) ), clinfo )
+ ELSEIF( PRESENT(pvdp2d) ) THEN
+ CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvdp2d(ix1:ix2,iy1:iy2) ), clinfo )
+ ELSEIF( PRESENT(pvsp3d) ) THEN
+ CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvsp3d(ix1:ix2,iy1:iy2,:) ), clinfo )
+ ELSEIF( PRESENT(pvdp3d) ) THEN
+ CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvdp3d(ix1:ix2,iy1:iy2,:) ), clinfo )
ENDIF
! add 1 to the size of the temporal dimension (not really useful...)
IF( iom_file(kiomid)%luld(idvar) ) iom_file(kiomid)%dimsz(iom_file(kiomid)%ndims(idvar), idvar) &
@@ -726,7 +754,7 @@ CONTAINS
IF(lwp) WRITE(numout,*) TRIM(clinfo)//' written ok'
ENDIF
!
- END SUBROUTINE iom_nf90_rp0123d_dp
+ END SUBROUTINE iom_nf90_rstput
SUBROUTINE iom_nf90_check( kstatus, cdinfo )
diff --git a/src/OCE/LBC/lib_mpp.F90 b/src/OCE/LBC/lib_mpp.F90
index 35d87d400327691fd366e0bc5e0932d8e9f73754..dfcb10a654ee4ded8ee82d80daa289b3187efef7 100644
--- a/src/OCE/LBC/lib_mpp.F90
+++ b/src/OCE/LBC/lib_mpp.F90
@@ -1234,7 +1234,7 @@ CONTAINS
+ aimag(ydda(ji)) + aimag(yddb(ji))
! The result is zt1 + zt2, after normalization.
- yddb(ji) = cmplx ( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1),wp )
+ yddb(ji) = cmplx ( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), dp )
END DO
!
END SUBROUTINE DDPDD_MPI
diff --git a/src/OCE/LDF/ldfslp.F90 b/src/OCE/LDF/ldfslp.F90
index 5a1fcf8ef324ef5d7cad1832d1df74cd7368b349..33d5957c26ad726ef3642871992dcd546c9ee425 100644
--- a/src/OCE/LDF/ldfslp.F90
+++ b/src/OCE/LDF/ldfslp.F90
@@ -610,7 +610,7 @@ CONTAINS
! END DO
! END DO
! END DO
-! CALL lbc_lnk( 'ldfslp', uslp , 'U', -1. ; CALL lbc_lnk( 'ldfslp', vslp , 'V', -1., wslpi, 'W', -1., wslpj, 'W', -1. )
+! CALL lbc_lnk( 'ldfslp', uslp , 'U', -1. ; CALL lbc_lnk( 'ldfslp', vslp , 'V', -1._wp, wslpi, 'W', -1._wp, wslpj, 'W', -1._wp )
!!gm ENDIF
ENDIF
!
diff --git a/src/OCE/SBC/sbcblk_algo_ecmwf.F90 b/src/OCE/SBC/sbcblk_algo_ecmwf.F90
index eb316c81a35d78dbbe837788e3d03a477cad8f03..252a19275feb0b2eba8ca4cb7cb3f5772b8cb9a3 100644
--- a/src/OCE/SBC/sbcblk_algo_ecmwf.F90
+++ b/src/OCE/SBC/sbcblk_algo_ecmwf.F90
@@ -417,27 +417,23 @@ CONTAINS
REAL(wp), DIMENSION(A2D(0)), INTENT(in) :: pzeta
!
INTEGER :: ji, jj ! dummy loop indices
- REAL(wp) :: zta, zx2, zx, ztmp, zpsi_unst, zpsi_stab, zstab, zc
+ REAL(wp) :: zta, zx2, zx, ztmp, zc
!!----------------------------------------------------------------------------------
zc = 5._wp/0.35_wp
DO_2D( 0, 0, 0, 0 )
!
zta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
-
- ! *** Unstable (Paulson 1970) [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
- zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 - 16z)^0.5
- zx = SQRT(zx2) ! (1 - 16z)^0.25
- ztmp = 1._wp + zx
- zpsi_unst = LOG( 0.125_wp*ztmp*ztmp*(1._wp + zx2) ) - 2._wp*ATAN( zx ) + 0.5_wp*rpi
-
- ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
- zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) &
- & - zta - 2._wp/3._wp*zc
!
- zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
- !
- psi_m_ecmwf(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable
- & + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable
+ IF (zta < 0.0_wp) THEN
+ ! *** Unstable (Paulson 1970) [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
+ zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 - 16z)^0.5
+ zx = SQRT(zx2) ! (1 - 16z)^0.25
+ ztmp = 1._wp + zx
+ psi_m_ecmwf(ji,jj) = LOG(0.125_wp*ztmp*ztmp*(1._wp + zx2)) - 2._wp*ATAN(zx) + 0.5_wp*rpi
+ ELSE
+ ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
+ psi_m_ecmwf(ji,jj) = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) - zta - 2._wp/3._wp*zc
+ ENDIF
!
END_2D
END FUNCTION psi_m_ecmwf
@@ -458,7 +454,7 @@ CONTAINS
REAL(wp), DIMENSION(A2D(0)), INTENT(in) :: pzeta
!
INTEGER :: ji, jj ! dummy loop indices
- REAL(wp) :: zta, zx2, zpsi_unst, zpsi_stab, zstab, zc
+ REAL(wp) :: zta, zx2, zc
!!----------------------------------------------------------------------------------
zc = 5._wp/0.35_wp
!
@@ -466,20 +462,15 @@ CONTAINS
!
zta = MIN(pzeta(ji,jj) , 5._wp) ! Very stable conditions (L positif and big!):
!
- ! *** Unstable (Paulson 1970) [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
- zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 -16z)^0.5
- zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
- !
- ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
- zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) &
+ IF (zta < 0.0_wp) THEN
+ ! *** Unstable (Paulson 1970) [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
+ zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 -16z)^0.5
+ psi_h_ecmwf(ji,jj) = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
+ ELSE
+ ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
+ psi_h_ecmwf(ji,jj) = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) &
& - ABS(1._wp + 2._wp/3._wp*zta)**1.5_wp - 2._wp/3._wp*zc + 1._wp
- !
- ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
- !
- zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
- !
- psi_h_ecmwf(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable
- & + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable
+ ENDIF
!
END_2D
END FUNCTION psi_h_ecmwf
diff --git a/src/OCE/SBC/sbcice_cice.F90 b/src/OCE/SBC/sbcice_cice.F90
index d05ccfd32cd33c0624cdb6633335f1dc7bbd4967..2af309172c40eed8d8cc9455f9510a1973e16781 100644
--- a/src/OCE/SBC/sbcice_cice.F90
+++ b/src/OCE/SBC/sbcice_cice.F90
@@ -175,7 +175,7 @@ CONTAINS
! there is no restart file.
! Values from a CICE restart file would overwrite this
IF( .NOT. ln_rstart ) THEN
- CALL nemo2cice( ts(:,:,1,jp_tem,Kmm) , sst , 'T' , 1.)
+ CALL nemo2cice( ts(:,:,1,jp_tem,Kmm) , sst , 'T' , 1._wp)
ENDIF
#endif
@@ -207,10 +207,10 @@ CONTAINS
fr_iu(:,:)=0.0
fr_iv(:,:)=0.0
- CALL cice2nemo(aice,fr_i, 'T', 1. )
+ CALL cice2nemo(aice,fr_i, 'T', 1._wp )
IF( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
DO jl=1,ncat
- CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
+ CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1._wp )
ENDDO
ENDIF
@@ -224,8 +224,8 @@ CONTAINS
CALL lbc_lnk( 'sbcice_cice', fr_iu , 'U', 1.0_wp, fr_iv , 'V', 1.0_wp )
! set the snow+ice mass
- CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
- CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
+ CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1._wp )
+ CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1._wp )
snwice_mass (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:) )
snwice_mass_b(:,:) = snwice_mass(:,:)
@@ -282,7 +282,7 @@ CONTAINS
ztmp(ji,jj) = 0.5 * ( fr_iu(ji,jj) * utau(ji,jj) &
+ fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
END_2D
- CALL nemo2cice(ztmp,strax,'F', -1. )
+ CALL nemo2cice(ztmp,strax,'F', -1._wp )
! y comp of wind stress (CI_2)
! V point to F point
@@ -290,7 +290,7 @@ CONTAINS
ztmp(ji,jj) = 0.5 * ( fr_iv(ji,jj) * vtau(ji,jj) &
+ fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
END_2D
- CALL nemo2cice(ztmp,stray,'F', -1. )
+ CALL nemo2cice(ztmp,stray,'F', -1._wp )
! Surface downward latent heat flux (CI_5)
IF(ksbc == jp_flx) THEN
@@ -316,7 +316,7 @@ CONTAINS
END_2D
ENDIF
DO jl=1,ncat
- CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
+ CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1._wp )
! GBM conductive flux through ice (CI_6)
! Convert to GBM
@@ -325,7 +325,7 @@ CONTAINS
ELSE
ztmp(:,:) = botmelt(:,:,jl)
ENDIF
- CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
+ CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1._wp )
! GBM surface heat flux (CI_7)
! Convert to GBM
@@ -334,7 +334,7 @@ CONTAINS
ELSE
ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
ENDIF
- CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
+ CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1._wp )
ENDDO
ELSE IF(ksbc == jp_blk) THEN
@@ -342,39 +342,39 @@ CONTAINS
! Pass bulk forcing fields to CICE (which will calculate heat fluxes etc itself)
! x comp and y comp of atmosphere surface wind (CICE expects on T points)
ztmp(:,:) = wndi_ice(:,:)
- CALL nemo2cice(ztmp,uatm,'T', -1. )
+ CALL nemo2cice(ztmp,uatm,'T', -1._wp )
ztmp(:,:) = wndj_ice(:,:)
- CALL nemo2cice(ztmp,vatm,'T', -1. )
+ CALL nemo2cice(ztmp,vatm,'T', -1._wp )
ztmp(:,:) = SQRT ( wndi_ice(:,:)**2 + wndj_ice(:,:)**2 )
- CALL nemo2cice(ztmp,wind,'T', 1. ) ! Wind speed (m/s)
+ CALL nemo2cice(ztmp,wind,'T', 1._wp ) ! Wind speed (m/s)
ztmp(:,:) = qsr_ice(:,:,1)
- CALL nemo2cice(ztmp,fsw,'T', 1. ) ! Incoming short-wave (W/m^2)
+ CALL nemo2cice(ztmp,fsw,'T', 1._wp ) ! Incoming short-wave (W/m^2)
ztmp(:,:) = qlw_ice(:,:,1)
- CALL nemo2cice(ztmp,flw,'T', 1. ) ! Incoming long-wave (W/m^2)
+ CALL nemo2cice(ztmp,flw,'T', 1._wp ) ! Incoming long-wave (W/m^2)
ztmp(:,:) = tatm_ice(:,:)
- CALL nemo2cice(ztmp,Tair,'T', 1. ) ! Air temperature (K)
- CALL nemo2cice(ztmp,potT,'T', 1. ) ! Potential temp (K)
+ CALL nemo2cice(ztmp,Tair,'T', 1._wp ) ! Air temperature (K)
+ CALL nemo2cice(ztmp,potT,'T', 1._wp ) ! Potential temp (K)
! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows
ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )
! Constant (101000.) atm pressure assumed
- CALL nemo2cice(ztmp,rhoa,'T', 1. ) ! Air density (kg/m^3)
+ CALL nemo2cice(ztmp,rhoa,'T', 1._wp ) ! Air density (kg/m^3)
ztmp(:,:) = qatm_ice(:,:)
- CALL nemo2cice(ztmp,Qa,'T', 1. ) ! Specific humidity (kg/kg)
+ CALL nemo2cice(ztmp,Qa,'T', 1._wp ) ! Specific humidity (kg/kg)
ztmp(:,:)=10.0
- CALL nemo2cice(ztmp,zlvl,'T', 1. ) ! Atmos level height (m)
+ CALL nemo2cice(ztmp,zlvl,'T', 1._wp ) ! Atmos level height (m)
! May want to check all values are physically realistic (as in CICE routine
! prepare_forcing)?
! Divide shortwave into spectral bands (as in prepare_forcing)
ztmp(:,:)=qsr_ice(:,:,1)*frcvdr ! visible direct
- CALL nemo2cice(ztmp,swvdr,'T', 1. )
+ CALL nemo2cice(ztmp,swvdr,'T', 1._wp )
ztmp(:,:)=qsr_ice(:,:,1)*frcvdf ! visible diffuse
- CALL nemo2cice(ztmp,swvdf,'T', 1. )
+ CALL nemo2cice(ztmp,swvdf,'T', 1._wp )
ztmp(:,:)=qsr_ice(:,:,1)*frcidr ! near IR direct
- CALL nemo2cice(ztmp,swidr,'T', 1. )
+ CALL nemo2cice(ztmp,swidr,'T', 1._wp )
ztmp(:,:)=qsr_ice(:,:,1)*frcidf ! near IR diffuse
- CALL nemo2cice(ztmp,swidf,'T', 1. )
+ CALL nemo2cice(ztmp,swidf,'T', 1._wp )
ENDIF
@@ -382,37 +382,37 @@ CONTAINS
! Ensure fsnow is positive (as in CICE routine prepare_forcing)
IF( iom_use('snowpre') ) CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit
ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)
- CALL nemo2cice(ztmp,fsnow,'T', 1. )
+ CALL nemo2cice(ztmp,fsnow,'T', 1._wp )
! Rainfall
IF( iom_use('precip') ) CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
- CALL nemo2cice(ztmp,frain,'T', 1. )
+ CALL nemo2cice(ztmp,frain,'T', 1._wp )
! Freezing/melting potential
! Calculated over NEMO leapfrog timestep (hence 2*dt)
nfrzmlt(:,:) = rho0 * rcp * e3t_m(:,:) * ( Tocnfrz-sst_m(:,:) ) / ( 2.0*dt )
ztmp(:,:) = nfrzmlt(:,:)
- CALL nemo2cice(ztmp,frzmlt,'T', 1. )
+ CALL nemo2cice(ztmp,frzmlt,'T', 1._wp )
! SST and SSS
- CALL nemo2cice(sst_m,sst,'T', 1. )
- CALL nemo2cice(sss_m,sss,'T', 1. )
+ CALL nemo2cice(sst_m,sst,'T', 1._wp )
+ CALL nemo2cice(sss_m,sss,'T', 1._wp )
! x comp and y comp of surface ocean current
! U point to F point
DO_2D( 1, 1, 1, 0 )
ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1)
END_2D
- CALL nemo2cice(ztmp,uocn,'F', -1. )
+ CALL nemo2cice(ztmp,uocn,'F', -1._wp )
! V point to F point
DO_2D( 1, 0, 1, 1 )
ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1)
END_2D
- CALL nemo2cice(ztmp,vocn,'F', -1. )
+ CALL nemo2cice(ztmp,vocn,'F', -1._wp )
IF( ln_ice_embd ) THEN !== embedded sea ice: compute representative ice top surface ==!
!
@@ -437,14 +437,14 @@ CONTAINS
ztmp(ji,jj)=0.5 * ( (zpice(ji+1,jj )-zpice(ji,jj )) * r1_e1u(ji,jj ) &
& + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1) ) * fmask(ji,jj,1)
END_2D
- CALL nemo2cice( ztmp,ss_tltx,'F', -1. )
+ CALL nemo2cice( ztmp,ss_tltx,'F', -1._wp )
! T point to F point
DO_2D( 1, 0, 1, 0 )
ztmp(ji,jj)=0.5 * ( (zpice(ji ,jj+1)-zpice(ji ,jj)) * r1_e2v(ji ,jj) &
& + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj) ) * fmask(ji,jj,1)
END_2D
- CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
+ CALL nemo2cice(ztmp,ss_tlty,'F', -1._wp )
!
END SUBROUTINE cice_sbc_in
@@ -466,7 +466,7 @@ CONTAINS
ENDIF
! x comp of ocean-ice stress
- CALL cice2nemo(strocnx,ztmp1,'F', -1. )
+ CALL cice2nemo(strocnx,ztmp1,'F', -1._wp )
ss_iou(:,:)=0.0
! F point to U point
DO_2D( 0, 0, 0, 0 )
@@ -475,7 +475,7 @@ CONTAINS
CALL lbc_lnk( 'sbcice_cice', ss_iou , 'U', -1.0_wp )
! y comp of ocean-ice stress
- CALL cice2nemo(strocny,ztmp1,'F', -1. )
+ CALL cice2nemo(strocny,ztmp1,'F', -1._wp )
ss_iov(:,:)=0.0
! F point to V point
@@ -494,8 +494,8 @@ CONTAINS
! Also need ice/ocean stress on T points so that taum can be updated
! This interpolation is already done in CICE so best to use those values
- CALL cice2nemo(strocnxT,ztmp1,'T',-1.)
- CALL cice2nemo(strocnyT,ztmp2,'T',-1.)
+ CALL cice2nemo(strocnxT,ztmp1,'T',-1._wp)
+ CALL cice2nemo(strocnyT,ztmp2,'T',-1._wp)
! Update taum with modulus of ice-ocean stress
! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here
@@ -518,11 +518,11 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
ENDIF
#if defined key_cice4
- CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
- CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
+ CALL cice2nemo(fresh_gbm,ztmp1,'T', 1._wp )
+ CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1._wp )
#else
- CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
- CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
+ CALL cice2nemo(fresh_ai,ztmp1,'T', 1._wp )
+ CALL cice2nemo(fsalt_ai,ztmp2,'T', 1._wp )
#endif
! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
@@ -556,9 +556,9 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
! Now add in ice / snow related terms
! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
#if defined key_cice4
- CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
+ CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1._wp )
#else
- CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
+ CALL cice2nemo(fswthru_ai,ztmp1,'T', 1._wp )
#endif
qsr(:,:)=qsr(:,:)+ztmp1(:,:)
CALL lbc_lnk( 'sbcice_cice', qsr , 'T', 1.0_wp )
@@ -568,9 +568,9 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
END_2D
#if defined key_cice4
- CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
+ CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1._wp )
#else
- CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
+ CALL cice2nemo(fhocn_ai,ztmp1,'T', 1._wp )
#endif
qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
@@ -578,10 +578,10 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
! Prepare for the following CICE time-step
- CALL cice2nemo(aice,fr_i,'T', 1. )
+ CALL cice2nemo(aice,fr_i,'T', 1._wp )
IF( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
DO jl=1,ncat
- CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
+ CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1._wp )
ENDDO
ENDIF
@@ -595,8 +595,8 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
CALL lbc_lnk( 'sbcice_cice', fr_iu , 'U', 1.0_wp, fr_iv , 'V', 1.0_wp )
! set the snow+ice mass
- CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
- CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
+ CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1._wp )
+ CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1._wp )
snwice_mass (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:) )
snwice_mass_b(:,:) = snwice_mass(:,:)
snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
@@ -628,16 +628,16 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
!
! x and y comp of ice velocity
!
- CALL cice2nemo(uvel,u_ice,'F', -1. )
- CALL cice2nemo(vvel,v_ice,'F', -1. )
+ CALL cice2nemo(uvel,u_ice,'F', -1._wp )
+ CALL cice2nemo(vvel,v_ice,'F', -1._wp )
!
! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out
!
! Snow and ice thicknesses (CO_2 and CO_3)
!
DO jl = 1, ncat
- CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1. )
- CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1. )
+ CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1._wp )
+ CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1._wp )
END DO
!
END SUBROUTINE cice_sbc_hadgam
diff --git a/src/OCE/TRA/traadv_fct.F90 b/src/OCE/TRA/traadv_fct.F90
index ce6121a6ebbf686a009ea5431b5c24b976d6d412..d9b3eb5858e865dacc8182506e46c6c2da77a28c 100644
--- a/src/OCE/TRA/traadv_fct.F90
+++ b/src/OCE/TRA/traadv_fct.F90
@@ -391,23 +391,25 @@ CONTAINS
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ikm1 ! local integer
- REAL(dp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars
- REAL(dp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
- REAL(dp), DIMENSION(T2D(nn_hls),jpk) :: zbetup, zbetdo, zbup, zbdo
+ REAL(wp) :: zpos, zneg, zbt, zbig ! local scalars
+ REAL(wp) :: zup, zdo ! - -
+ REAL(wp), DIMENSION(T2D(nn_hls),jpk) :: zbetup, zbetdo, zbup, zbdo
!!----------------------------------------------------------------------
!
- zbig = 1.e+40_dp
- zrtrn = 1.e-15_dp
- zbetup(:,:,:) = 0._dp ; zbetdo(:,:,:) = 0._dp
+ zbig = HUGE(1._wp)
+ zbetup(:,:,jpk) = zbig ; zbetdo(:,:,jpk) = zbig
! Search local extrema
! --------------------
! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
- zbup(ji,jj,jk) = MAX( pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ), &
- & paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ) )
- zbdo(ji,jj,jk) = MIN( pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ), &
- & paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ) )
+ IF( tmask(ji,jj,jk) == 1._wp ) THEN
+ zbup(ji,jj,jk) = MAX( pbef(ji,jj,jk), paft(ji,jj,jk) )
+ zbdo(ji,jj,jk) = MIN( pbef(ji,jj,jk), paft(ji,jj,jk) )
+ ELSE
+ zbup(ji,jj,jk) = -zbig
+ zbdo(ji,jj,jk) = zbig
+ ENDIF
END_3D
DO jk = 1, jpkm1
@@ -427,41 +429,38 @@ CONTAINS
& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
! positive part of the flux
- zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
- & + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
- & + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
+ zpos = MAX( 0._wp, paa(ji-1,jj ,jk ) ) - MIN( 0._wp, paa(ji ,jj ,jk ) ) &
+ & + MAX( 0._wp, pbb(ji ,jj-1,jk ) ) - MIN( 0._wp, pbb(ji ,jj ,jk ) ) &
+ & + MAX( 0._wp, pcc(ji ,jj ,jk+1) ) - MIN( 0._wp, pcc(ji ,jj ,jk ) )
! negative part of the flux
- zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
- & + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
- & + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
+ zneg = MAX( 0._wp, paa(ji ,jj ,jk ) ) - MIN( 0._wp, paa(ji-1,jj ,jk ) ) &
+ & + MAX( 0._wp, pbb(ji ,jj ,jk ) ) - MIN( 0._wp, pbb(ji ,jj-1,jk ) ) &
+ & + MAX( 0._wp, pcc(ji ,jj ,jk ) ) - MIN( 0._wp, pcc(ji ,jj ,jk+1) )
! up & down beta terms
zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kaa) / p2dt
- zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
- zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
+ IF( zup /= -zbig .AND. zpos /= 0._wp ) THEN ; zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / zpos * zbt
+ ELSE ; zbetup(ji,jj,jk) = zbig
+ ENDIF
+ IF( zdo /= zbig .AND. zneg /= 0._wp ) THEN ; zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / zneg * zbt
+ ELSE ; zbetdo(ji,jj,jk) = zbig
+ ENDIF
END_2D
END DO
- ! 3. monotonic flux in the i & j direction (paa & pbb)
- ! ----------------------------------------
+ ! 3. monotonic flux in the i, j and k direction (paa, pbb and pcc)
+ ! ----------------------------------------------------------------
DO_3D( 1, 0, 1, 0, 1, jpkm1 )
- zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
- zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
- zcu = ( 0.5 + SIGN( 0.5_wp , paa(ji,jj,jk) ) )
- paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
-
- zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
- zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
- zcv = ( 0.5 + SIGN( 0.5_wp , pbb(ji,jj,jk) ) )
- pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
-
- ! monotonic flux in the k direction, i.e. pcc
- ! -------------------------------------------
- za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
- zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
- zc = ( 0.5 + SIGN( 0.5_wp , pcc(ji,jj,jk+1) ) )
- pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
+ IF( paa(ji,jj,jk) > 0._wp ) THEN ; paa(ji,jj,jk) = paa(ji,jj,jk) * MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
+ ELSE ; paa(ji,jj,jk) = paa(ji,jj,jk) * MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
+ ENDIF
+ IF( pbb(ji,jj,jk) > 0._wp ) THEN ; pbb(ji,jj,jk) = pbb(ji,jj,jk) * MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
+ ELSE ; pbb(ji,jj,jk) = pbb(ji,jj,jk) * MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
+ ENDIF
+ IF( pcc(ji,jj,jk+1) > 0._wp ) THEN ; pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * MIN( 1._wp, zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
+ ELSE ; pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * MIN( 1._wp, zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
+ ENDIF
END_3D
!
END SUBROUTINE nonosc
diff --git a/src/OCE/TRA/traadv_ubs.F90 b/src/OCE/TRA/traadv_ubs.F90
index 02d4d064912b58a344d7b7a1e135b7d66eca3bef..0ebea711f769f6bc2e07a42986f7089ed3a021e6 100644
--- a/src/OCE/TRA/traadv_ubs.F90
+++ b/src/OCE/TRA/traadv_ubs.F90
@@ -285,7 +285,7 @@ CONTAINS
REAL(wp), DIMENSION(T2D(nn_hls),jpk) :: zbetup, zbetdo ! 3D workspace
!!----------------------------------------------------------------------
!
- zbig = 1.e+38_wp
+ zbig = 1.e+20_wp ! works ok with simple/double precison
zrtrn = 1.e-15_wp
zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp
!
diff --git a/src/OCE/TRA/traadv_ubs_lf.F90 b/src/OCE/TRA/traadv_ubs_lf.F90
index 4a2b696c755c35e9ca862ab4c7eb6522102a6d9b..56b5e26b6ede932b6dc51a1a01c0b8bcc1b1c364 100644
--- a/src/OCE/TRA/traadv_ubs_lf.F90
+++ b/src/OCE/TRA/traadv_ubs_lf.F90
@@ -298,7 +298,7 @@ CONTAINS
REAL(wp), DIMENSION(T2D(nn_hls),jpk) :: zbetup, zbetdo ! 3D workspace
!!----------------------------------------------------------------------
!
- zbig = 1.e+38_wp
+ zbig = 1.e+20_wp ! works ok with simple/double precison
zrtrn = 1.e-15_wp
zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp
!
diff --git a/src/OCE/TRA/traqsr.F90 b/src/OCE/TRA/traqsr.F90
index 765d7313e72ffba4c2e9facd34b42e58e1af38c0..0975038a9980457282ac657899511c3d80274c86 100644
--- a/src/OCE/TRA/traqsr.F90
+++ b/src/OCE/TRA/traqsr.F90
@@ -77,6 +77,7 @@ MODULE traqsr
REAL(wp) :: r1_si0 ! all schemes : infrared = 1/rn_si0
REAL(wp) :: r1_si1 ! 2 band : mean RGB = 1/rn_si1
REAL(wp) :: r1_LR, r1_LG, r1_LB ! RGB with constant Chl (np_RGB)
+ REAL(wp) :: zz0
!
REAL(wp) , PUBLIC, DIMENSION(3,61) :: rkrgb ! tabulated attenuation coefficients for RGB absorption
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read)
@@ -209,8 +210,9 @@ CONTAINS
!
! sea-ice: store the 1st ocean level attenuation coefficient
DO_2D( 0, 0, 0, 0 )
- IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
- ELSE ; fraqsr_1lev(ji,jj) = 1._wp
+ zz0 = r1_rho0_rcp * qsr(ji,jj) ! test zz0 and not qsr for rounding errors in single precision
+ IF( zz0 /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / zz0
+ ELSE ; fraqsr_1lev(ji,jj) = 1._wp
ENDIF
END_2D
!
diff --git a/src/OCE/ZDF/zdfddm.F90 b/src/OCE/ZDF/zdfddm.F90
index 5fa6c719accb7b97fe379cb548384b3245525110..6d02b23ce600b82f77031c6fbf6f3ea6c76b4463 100644
--- a/src/OCE/ZDF/zdfddm.F90
+++ b/src/OCE/ZDF/zdfddm.F90
@@ -142,7 +142,7 @@ CONTAINS
! salt fingering
zrr = zrau(ji,jj) / rn_hsbfr
zrr = zrr * zrr
- zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
+ zavfs = rn_avts / ( 1. + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
zavft = 0.7 * zavfs * zinr
! diffusive layering
zavdt = 1.3635e-6 * EXP( 4.6 * EXP( -0.54*(zinr-1.) ) ) * zmsks(ji,jj) * zmskd1(ji,jj)
diff --git a/src/OCE/ZDF/zdftke.F90 b/src/OCE/ZDF/zdftke.F90
index 4435eb13ee1157b1542b58e7bf9444ed9080a3bc..f2499f8bce6d3f9efab346bfefe772db2ee2baca 100644
--- a/src/OCE/ZDF/zdftke.F90
+++ b/src/OCE/ZDF/zdftke.F90
@@ -218,6 +218,7 @@ CONTAINS
REAL(wp), DIMENSION(T2D(0)) :: zice_fra, zhlc, zus3, zWlc2
REAL(wp), DIMENSION(T2D(0),jpk) :: zpelc, zdiag, zd_up, zd_lw
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztmp ! for diags
+ REAL(wp) :: zdiv
!!--------------------------------------------------------------------
!
zbbrau = rn_ebb / rho0 ! Local constant initialisation
@@ -369,7 +370,16 @@ CONTAINS
IF (rn2b(ji,jj,jk) <= 0.0_wp) then
zri = 0.0_wp
ELSE
- zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
+ ! This logic is to avoid divide-by-zero errors which can occur for single-precision
+ ! The actual value you choose for the denominator instead of zero doesn't really
+ ! matter, as long as it is very small and so triggers the same logic below with the
+ ! inverse Prandtl number
+ zdiv = p_sh2(ji,jj,jk) + rn_bshear
+ IF (zdiv == 0.0_wp) THEN
+ zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / rn_bshear
+ ELSE
+ zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / zdiv
+ ENDIF
ENDIF
! ! inverse of Prandtl number
apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) )
diff --git a/src/OCE/lib_fortran.F90 b/src/OCE/lib_fortran.F90
index e7538e7c89cd3ee0a45db16bda946154de5e0ad3..0a4ac708ded389bce46657c7d4eb9ceef643c8c4 100644
--- a/src/OCE/lib_fortran.F90
+++ b/src/OCE/lib_fortran.F90
@@ -241,7 +241,7 @@ CONTAINS
& + AIMAG(ydda) + AIMAG(yddb)
!
! The result is t1 + t2, after normalization.
- yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), wp )
+ yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), dp )
!
END SUBROUTINE DDPDD
diff --git a/src/OCE/module_example.F90 b/src/OCE/module_example.F90
index b0198f02333a07c85e03c7f7b54a742bca946536..fc9c8c8c9b334a5f1c6abefeb9e9356f6996694a 100644
--- a/src/OCE/module_example.F90
+++ b/src/OCE/module_example.F90
@@ -126,7 +126,7 @@ CONTAINS
!
! WARNING! the lbc_lnk call could not be compatible with the tiling approach
! please refer to the manual for how to adapt your code
- CALL lbc_lnk( 'module_example', avm, 'T', 1. ) ! Lateral boundary conditions (unchanged sign)
+ CALL lbc_lnk( 'module_example', avm, 'T', 1._wp ) ! Lateral boundary conditions (unchanged sign)
!
END SUBROUTINE exa_mpl
diff --git a/src/OCE/stpmlf.F90 b/src/OCE/stpmlf.F90
index 6a49213c47f818a2189ac110fa366df8d017f667..9d7c7c0adff311c450dc5c9bdccaa8cde095c9c9 100644
--- a/src/OCE/stpmlf.F90
+++ b/src/OCE/stpmlf.F90
@@ -538,8 +538,8 @@ CONTAINS
CALL Agrif_dyn( kt )
# endif
! ! local domain boundaries (T-point, unchanged sign)
- CALL lbc_lnk( 'finalize_lbc', puu(:,:,:, Kaa), 'U', -1., pvv(:,:,: ,Kaa), 'V', -1. &
- & , pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )
+ CALL lbc_lnk( 'finalize_lbc', puu(:,:,:, Kaa), 'U', -1._wp, pvv(:,:,: ,Kaa), 'V', -1._wp &
+ & , pts(:,:,:,jp_tem,Kaa), 'T', 1._wp, pts(:,:,:,jp_sal,Kaa), 'T', 1._wp )
!
! lbc_lnk needed for zdf_sh2, moved here to allow tiling in zdf_phy
IF( l_zdfsh2 ) CALL lbc_lnk( 'stp', avm_k, 'W', 1.0_wp )
diff --git a/src/OCE/stprk3_stg.F90 b/src/OCE/stprk3_stg.F90
index d1c35aa965a955da0346005d378d7f854e9ae0db..d45a6c744a607faa3fd2312de00c2ccf694e299d 100644
--- a/src/OCE/stprk3_stg.F90
+++ b/src/OCE/stprk3_stg.F90
@@ -447,8 +447,8 @@ CONTAINS
CALL Agrif_dyn( kstp, kstg )
# endif
! !* local domain boundaries (T-point, unchanged sign)
- CALL lbc_lnk( 'stp_RK3_stg', uu(:,:,:, Kaa), 'U', -1., vv(:,:,: ,Kaa), 'V', -1. &
- & , ts(:,:,:,jp_tem,Kaa), 'T', 1., ts(:,:,:,jp_sal,Kaa), 'T', 1. )
+ CALL lbc_lnk( 'stp_RK3_stg', uu(:,:,:, Kaa), 'U', -1._wp, vv(:,:,: ,Kaa), 'V', -1._wp &
+ & , ts(:,:,:,jp_tem,Kaa), 'T', 1._wp, ts(:,:,:,jp_sal,Kaa), 'T', 1._wp )
!
IF( l_zdfsh2 ) CALL lbc_lnk( 'stp_RK3_stg', avm_k, 'W', 1.0_wp ) ! lbc_lnk needed for zdf_sh2, moved here to allow tiling in zdf_phy
!
diff --git a/src/OCE/timing.F90 b/src/OCE/timing.F90
index c63003d8b66e5ffb901c471502339be6d98c2eeb..4c1becfe14fdbb9d97515b9e884015d8842ccfa4 100644
--- a/src/OCE/timing.F90
+++ b/src/OCE/timing.F90
@@ -39,7 +39,7 @@ MODULE timing
CHARACTER(LEN=20) :: cname
CHARACTER(LEN=20) :: surname
INTEGER :: rank
- REAL(wp) :: t_cpu, t_clock, tsum_cpu, tsum_clock, tmax_cpu, tmax_clock, tmin_cpu, tmin_clock, tsub_cpu, tsub_clock
+ REAL(dp) :: t_cpu, t_clock, tsum_cpu, tsum_clock, tmax_cpu, tmax_clock, tmin_cpu, tmin_clock, tsub_cpu, tsub_clock
INTEGER :: ncount, ncount_max, ncount_rate
INTEGER :: niter
LOGICAL :: l_tdone
@@ -50,8 +50,8 @@ MODULE timing
TYPE alltimer
CHARACTER(LEN=20), DIMENSION(:), POINTER :: cname => NULL()
- REAL(wp), DIMENSION(:), POINTER :: tsum_cpu => NULL()
- REAL(wp), DIMENSION(:), POINTER :: tsum_clock => NULL()
+ REAL(dp), DIMENSION(:), POINTER :: tsum_cpu => NULL()
+ REAL(dp), DIMENSION(:), POINTER :: tsum_clock => NULL()
INTEGER, DIMENSION(:), POINTER :: niter => NULL()
TYPE(alltimer), POINTER :: next => NULL()
TYPE(alltimer), POINTER :: prev => NULL()
@@ -62,14 +62,14 @@ MODULE timing
TYPE(timer), POINTER :: s_timer_old => NULL()
TYPE(timer), POINTER :: s_wrk => NULL()
- REAL(wp) :: t_overclock, t_overcpu
+ REAL(dp) :: t_overclock, t_overcpu
LOGICAL :: l_initdone = .FALSE.
INTEGER :: nsize
! Variables for coarse grain timing
- REAL(wp) :: tot_etime, tot_ctime
- REAL(kind=wp), DIMENSION(2) :: t_elaps, t_cpu
- REAL(wp), ALLOCATABLE, DIMENSION(:) :: all_etime, all_ctime
+ REAL(dp) :: tot_etime, tot_ctime
+ REAL(kind=dp), DIMENSION(2) :: t_elaps, t_cpu
+ REAL(dp), ALLOCATABLE, DIMENSION(:) :: all_etime, all_ctime
INTEGER :: nfinal_count, ncount, ncount_rate, ncount_max
INTEGER, DIMENSION(8) :: nvalues
CHARACTER(LEN=8), DIMENSION(2) :: cdate
@@ -136,7 +136,7 @@ CONTAINS
CHARACTER(len=*), INTENT(in), OPTIONAL :: csection
!
INTEGER :: ifinal_count, iperiods
- REAL(wp) :: zcpu_end, zmpitime,zcpu_raw,zclock_raw
+ REAL(dp) :: zcpu_end, zmpitime,zcpu_raw,zclock_raw
!
s_wrk => NULL()
@@ -218,7 +218,7 @@ CONTAINS
!! ** Purpose : open timing output file
!!----------------------------------------------------------------------
INTEGER :: iperiods, istart_count, ifinal_count
- REAL(wp) :: zdum
+ REAL(dp) :: zdum
LOGICAL :: ll_f
CHARACTER(len=*), INTENT(in), OPTIONAL :: clname
CHARACTER(len=20) :: cln
@@ -295,9 +295,9 @@ CONTAINS
INTEGER :: ji
LOGICAL :: ll_ord, ll_averep
CHARACTER(len=120) :: clfmt
- REAL(wp), DIMENSION(:), ALLOCATABLE :: timing_glob
- REAL(wp) :: zsypd ! simulated years per day (Balaji 2017)
- REAL(wp) :: zperc, ztot
+ REAL(dp), DIMENSION(:), ALLOCATABLE :: timing_glob
+ REAL(dp) :: zsypd ! simulated years per day (Balaji 2017)
+ REAL(dp) :: zperc, ztot
ll_averep = .TRUE.
@@ -387,16 +387,16 @@ CONTAINS
WRITE(numtime,*) ' warning: includes restarts writing time if output before nitend... '
WRITE(numtime,*) ' '
DO ji = 1, jpnij
- zperc = 0._wp ; zsypd = 0._wp
+ zperc = 0._dp ; zsypd = 0._dp
ztot = SUM( timing_glob(4*ji-3:4*ji-1) )
WRITE(numtime,'(A28,F11.6, A34,I8)') 'Computing time : ',timing_glob(4*ji-3), ' on MPI rank : ', ji
- IF ( ztot /= 0._wp ) zperc = timing_glob(4*ji-2) / ztot * 100.
+ IF ( ztot /= 0._dp ) zperc = timing_glob(4*ji-2) / ztot * 100.
WRITE(numtime,'(A28,F11.6,A2, F4.1,A3,A25,I8)') 'Waiting lbc_lnk time : ',timing_glob(4*ji-2) &
& , ' (', zperc,' %)', ' on MPI rank : ', ji
- IF ( ztot /= 0._wp ) zperc = timing_glob(4*ji-1) / ztot * 100.
+ IF ( ztot /= 0._dp ) zperc = timing_glob(4*ji-1) / ztot * 100.
WRITE(numtime,'(A28,F11.6,A2, F4.1,A3,A25,I8)') 'Waiting global time : ',timing_glob(4*ji-1) &
& , ' (', zperc,' %)', ' on MPI rank : ', ji
- IF ( timing_glob(4*ji) /= 0._wp ) zsypd = rn_Dt * REAL(nitend-nit000-1, wp) / (timing_glob(4*ji) * 365.)
+ IF ( timing_glob(4*ji) /= 0._dp ) zsypd = rn_Dt * REAL(nitend-nit000-1, dp) / (timing_glob(4*ji) * 365.)
WRITE(numtime,'(A28,F11.6,A7,F10.3,A2,A15,I8)') 'Total time : ',timing_glob(4*ji ) &
& , ' (SYPD: ', zsypd, ')', ' on MPI rank : ', ji
END DO
@@ -447,7 +447,7 @@ CONTAINS
s_timer => s_timer_root
clfmt = '(1x,a,4x,f12.3,6x,f12.3,x,f12.3,2x,f12.3,6x,f7.3,2x,i9)'
DO WHILE ( ASSOCIATED(s_timer) )
- IF( s_timer%tsum_clock > 0._wp ) &
+ IF( s_timer%tsum_clock > 0._dp ) &
WRITE(numtime,TRIM(clfmt)) s_timer%cname, &
& s_timer%tsum_clock,s_timer%tsum_clock*100./t_elaps(2), &
& s_timer%tsum_cpu ,s_timer%tsum_cpu*100./t_cpu(2) , &
@@ -489,8 +489,8 @@ CONTAINS
RETURN
END IF
sl_timer_glob_root%cname(:) = ''
- sl_timer_glob_root%tsum_cpu(:) = 0._wp
- sl_timer_glob_root%tsum_clock(:) = 0._wp
+ sl_timer_glob_root%tsum_cpu(:) = 0._dp
+ sl_timer_glob_root%tsum_clock(:) = 0._dp
sl_timer_glob_root%niter(:) = 0
sl_timer_glob_root%next => NULL()
sl_timer_glob_root%prev => NULL()
@@ -506,16 +506,16 @@ CONTAINS
IF( narea .EQ. 1 ) THEN
ALLOCATE(sl_timer_ave_root)
sl_timer_ave_root%cname = ''
- sl_timer_ave_root%t_cpu = 0._wp
- sl_timer_ave_root%t_clock = 0._wp
- sl_timer_ave_root%tsum_cpu = 0._wp
- sl_timer_ave_root%tsum_clock = 0._wp
- sl_timer_ave_root%tmax_cpu = 0._wp
- sl_timer_ave_root%tmax_clock = 0._wp
- sl_timer_ave_root%tmin_cpu = 0._wp
- sl_timer_ave_root%tmin_clock = 0._wp
- sl_timer_ave_root%tsub_cpu = 0._wp
- sl_timer_ave_root%tsub_clock = 0._wp
+ sl_timer_ave_root%t_cpu = 0._dp
+ sl_timer_ave_root%t_clock = 0._dp
+ sl_timer_ave_root%tsum_cpu = 0._dp
+ sl_timer_ave_root%tsum_clock = 0._dp
+ sl_timer_ave_root%tmax_cpu = 0._dp
+ sl_timer_ave_root%tmax_clock = 0._dp
+ sl_timer_ave_root%tmin_cpu = 0._dp
+ sl_timer_ave_root%tmin_clock = 0._dp
+ sl_timer_ave_root%tsub_cpu = 0._dp
+ sl_timer_ave_root%tsub_clock = 0._dp
sl_timer_ave_root%ncount = 0
sl_timer_ave_root%ncount_rate = 0
sl_timer_ave_root%ncount_max = 0
@@ -642,10 +642,10 @@ CONTAINS
!
INTEGER :: idum, icode
INTEGER, ALLOCATABLE, DIMENSION(:) :: iall_rank
- REAL(wp) :: ztot_ratio
- REAL(wp) :: zmax_etime, zmax_ctime, zmax_ratio, zmin_etime, zmin_ctime, zmin_ratio
- REAL(wp) :: zavg_etime, zavg_ctime, zavg_ratio
- REAL(wp), ALLOCATABLE, DIMENSION(:) :: zall_ratio
+ REAL(dp) :: ztot_ratio
+ REAL(dp) :: zmax_etime, zmax_ctime, zmax_ratio, zmin_etime, zmin_ctime, zmin_ratio
+ REAL(dp) :: zavg_etime, zavg_ctime, zavg_ratio
+ REAL(dp), ALLOCATABLE, DIMENSION(:) :: zall_ratio
CHARACTER(LEN=128), dimension(8) :: cllignes
CHARACTER(LEN=128) :: clhline, clstart_date, clfinal_date
CHARACTER(LEN=2048) :: clfmt
@@ -656,19 +656,19 @@ CONTAINS
iall_rank(:) = (/ (idum,idum=0,jpnij-1) /)
! Compute elapse user time
- zavg_etime = tot_etime/REAL(jpnij,wp)
+ zavg_etime = tot_etime/REAL(jpnij,dp)
zmax_etime = MAXVAL(all_etime(:))
zmin_etime = MINVAL(all_etime(:))
! Compute CPU user time
- zavg_ctime = tot_ctime/REAL(jpnij,wp)
+ zavg_ctime = tot_ctime/REAL(jpnij,dp)
zmax_ctime = MAXVAL(all_ctime(:))
zmin_ctime = MINVAL(all_ctime(:))
! Compute cpu/elapsed ratio
zall_ratio(:) = all_ctime(:) / all_etime(:)
ztot_ratio = SUM(all_ctime(:))/SUM(all_etime(:))
- zavg_ratio = SUM(zall_ratio(:))/REAL(jpnij,wp)
+ zavg_ratio = SUM(zall_ratio(:))/REAL(jpnij,dp)
zmax_ratio = MAXVAL(zall_ratio(:))
zmin_ratio = MINVAL(zall_ratio(:))
@@ -714,16 +714,16 @@ CONTAINS
IF( .NOT. ASSOCIATED(s_timer_root) ) THEN
ALLOCATE(s_timer_root)
s_timer_root%cname = cdinfo
- s_timer_root%t_cpu = 0._wp
- s_timer_root%t_clock = 0._wp
- s_timer_root%tsum_cpu = 0._wp
- s_timer_root%tsum_clock = 0._wp
- s_timer_root%tmax_cpu = 0._wp
- s_timer_root%tmax_clock = 0._wp
- s_timer_root%tmin_cpu = 0._wp
- s_timer_root%tmin_clock = 0._wp
- s_timer_root%tsub_cpu = 0._wp
- s_timer_root%tsub_clock = 0._wp
+ s_timer_root%t_cpu = 0._dp
+ s_timer_root%t_clock = 0._dp
+ s_timer_root%tsum_cpu = 0._dp
+ s_timer_root%tsum_clock = 0._dp
+ s_timer_root%tmax_cpu = 0._dp
+ s_timer_root%tmax_clock = 0._dp
+ s_timer_root%tmin_cpu = 0._dp
+ s_timer_root%tmin_clock = 0._dp
+ s_timer_root%tsub_cpu = 0._dp
+ s_timer_root%tsub_clock = 0._dp
s_timer_root%ncount = 0
s_timer_root%ncount_rate = 0
s_timer_root%ncount_max = 0
@@ -738,16 +738,16 @@ CONTAINS
!
ALLOCATE(s_timer_old)
s_timer_old%cname = cdinfo
- s_timer_old%t_cpu = 0._wp
- s_timer_old%t_clock = 0._wp
- s_timer_old%tsum_cpu = 0._wp
- s_timer_old%tsum_clock = 0._wp
- s_timer_old%tmax_cpu = 0._wp
- s_timer_old%tmax_clock = 0._wp
- s_timer_old%tmin_cpu = 0._wp
- s_timer_old%tmin_clock = 0._wp
- s_timer_old%tsub_cpu = 0._wp
- s_timer_old%tsub_clock = 0._wp
+ s_timer_old%t_cpu = 0._dp
+ s_timer_old%t_clock = 0._dp
+ s_timer_old%tsum_cpu = 0._dp
+ s_timer_old%tsum_clock = 0._dp
+ s_timer_old%tmax_cpu = 0._dp
+ s_timer_old%tmax_clock = 0._dp
+ s_timer_old%tmin_cpu = 0._dp
+ s_timer_old%tmin_clock = 0._dp
+ s_timer_old%tsub_cpu = 0._dp
+ s_timer_old%tsub_clock = 0._dp
s_timer_old%ncount = 0
s_timer_old%ncount_rate = 0
s_timer_old%ncount_max = 0
@@ -780,16 +780,16 @@ CONTAINS
! write(*,*) 'after allocation of next'
s_timer%next%cname = cdinfo
- s_timer%next%t_cpu = 0._wp
- s_timer%next%t_clock = 0._wp
- s_timer%next%tsum_cpu = 0._wp
- s_timer%next%tsum_clock = 0._wp
- s_timer%next%tmax_cpu = 0._wp
- s_timer%next%tmax_clock = 0._wp
- s_timer%next%tmin_cpu = 0._wp
- s_timer%next%tmin_clock = 0._wp
- s_timer%next%tsub_cpu = 0._wp
- s_timer%next%tsub_clock = 0._wp
+ s_timer%next%t_cpu = 0._dp
+ s_timer%next%t_clock = 0._dp
+ s_timer%next%tsum_cpu = 0._dp
+ s_timer%next%tsum_clock = 0._dp
+ s_timer%next%tmax_cpu = 0._dp
+ s_timer%next%tmax_clock = 0._dp
+ s_timer%next%tmin_cpu = 0._dp
+ s_timer%next%tmin_clock = 0._dp
+ s_timer%next%tsub_cpu = 0._dp
+ s_timer%next%tsub_clock = 0._dp
s_timer%next%ncount = 0
s_timer%next%ncount_rate = 0
s_timer%next%ncount_max = 0
diff --git a/src/SWE/stpmlf.F90 b/src/SWE/stpmlf.F90
index 4c542fb18111582d1888d1f036758de17b59b3ea..f358a28269a2aa93c1a6b199071263b839c606ef 100644
--- a/src/SWE/stpmlf.F90
+++ b/src/SWE/stpmlf.F90
@@ -196,9 +196,9 @@ CONTAINS
ENDIF
ENDIF
- CALL lbc_lnk( 'stp_MLF', uu (:,:,:,Nnn), 'U', -1., vv (:,:,:,Nnn), 'V', -1., & !* local domain boundaries
- & uu (:,:,:,Naa), 'U', -1., vv (:,:,:,Naa), 'V', -1. )
- CALL lbc_lnk( 'stp_MLF', r3u(:,:, Naa), 'U', 1., r3v(:,:, Naa), 'V', 1. )
+ CALL lbc_lnk( 'stp_MLF', uu (:,:,:,Nnn), 'U', -1._wp, vv (:,:,:,Nnn), 'V', -1._wp, & !* local domain boundaries
+ & uu (:,:,:,Naa), 'U', -1._wp, vv (:,:,:,Naa), 'V', -1._wp )
+ CALL lbc_lnk( 'stp_MLF', r3u(:,:, Naa), 'U', 1._wp, r3v(:,:, Naa), 'V', 1._wp )
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Set boundary conditions, time filter and swap time levels
diff --git a/src/SWE/stprk3.F90 b/src/SWE/stprk3.F90
index 0aa6cd2f7a22b14fa7ddf69adb2abaec2a5d975a..d0597b4451d55f0d433a4be1e95da292ddd00294 100644
--- a/src/SWE/stprk3.F90
+++ b/src/SWE/stprk3.F90
@@ -170,8 +170,8 @@ CONTAINS
END_3D
ENDIF
!
- CALL lbc_lnk( 'stp_RK3', uu (:,:,:,Naa), 'U', -1., vv (:,:,:,Naa), 'V', -1. )
- CALL lbc_lnk( 'stp_RK3', r3u(:,:, Naa), 'U', 1., r3v(:,:, Naa), 'V', 1. )
+ CALL lbc_lnk( 'stp_RK3', uu (:,:,:,Naa), 'U', -1._wp, vv (:,:,:,Naa), 'V', -1._wp )
+ CALL lbc_lnk( 'stp_RK3', r3u(:,:, Naa), 'U', 1._wp, r3v(:,:, Naa), 'V', 1._wp )
!
! !== Swap time levels ==!
Nrhs= Nnn
@@ -236,8 +236,8 @@ CONTAINS
END_3D
ENDIF
!
- CALL lbc_lnk( 'stp_RK3', uu (:,:,:,Naa), 'U', -1., vv (:,:,:,Naa), 'V', -1. )
- CALL lbc_lnk( 'stp_RK3', r3u(:,:, Naa), 'U', 1., r3v(:,:, Naa), 'V', 1. )
+ CALL lbc_lnk( 'stp_RK3', uu (:,:,:,Naa), 'U', -1._wp, vv (:,:,:,Naa), 'V', -1._wp )
+ CALL lbc_lnk( 'stp_RK3', r3u(:,:, Naa), 'U', 1._wp, r3v(:,:, Naa), 'V', 1._wp )
!
! !== Swap time levels ==!
Nrhs= Nnn
@@ -300,8 +300,8 @@ CONTAINS
END_3D
ENDIF
!
- CALL lbc_lnk( 'stp_RK3', uu (:,:,:,Naa), 'U', -1., vv (:,:,:,Naa), 'V', -1. )
- CALL lbc_lnk( 'stp_RK3', r3u(:,:, Naa), 'U', 1., r3v(:,:, Naa), 'V', 1. )
+ CALL lbc_lnk( 'stp_RK3', uu (:,:,:,Naa), 'U', -1._wp, vv (:,:,:,Naa), 'V', -1._wp )
+ CALL lbc_lnk( 'stp_RK3', r3u(:,:, Naa), 'U', 1._wp, r3v(:,:, Naa), 'V', 1._wp )
!
! !== Swap time levels ==!
!
diff --git a/src/TOP/TRP/trdmxl_trc.F90 b/src/TOP/TRP/trdmxl_trc.F90
index 1d13a2e85c82820a74218d60de645400a3a02f6f..80f20b0fbbbf54a1a97d1998aed01e625ce3fed8 100644
--- a/src/TOP/TRP/trdmxl_trc.F90
+++ b/src/TOP/TRP/trdmxl_trc.F90
@@ -289,7 +289,7 @@ CONTAINS
DO jn = 1, jptra
IF( ln_trdtrc(jn) ) THEN
DO jl = 1, jpltrd_trc
- CALL lbc_lnk( 'trdmxl_trc', tmltrd_trc(:,:,jl,jn), 'T', 1. ) ! lateral boundary conditions
+ CALL lbc_lnk( 'trdmxl_trc', tmltrd_trc(:,:,jl,jn), 'T', 1._wp ) ! lateral boundary conditions
END DO
ENDIF
END DO
@@ -424,8 +424,8 @@ CONTAINS
!-- Lateral boundary conditions
IF ( cn_cfg .NE. 'gyre' ) THEN
- CALL lbc_lnk( 'trdmxl_trc', ztmltot(:,:,jn) , 'T', 1. , ztmlres(:,:,jn) , 'T', 1., &
- & ztmlatf(:,:,jn) , 'T', 1. , ztmlrad(:,:,jn) , 'T', 1. )
+ CALL lbc_lnk( 'trdmxl_trc', ztmltot(:,:,jn) , 'T', 1._wp , ztmlres(:,:,jn) , 'T', 1._wp, &
+ & ztmlatf(:,:,jn) , 'T', 1._wp , ztmlrad(:,:,jn) , 'T', 1._wp )
ENDIF
@@ -475,9 +475,9 @@ CONTAINS
!-- Lateral boundary conditions
IF ( cn_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
- CALL lbc_lnk( 'trdmxl_trc', ztmltot2(:,:,jn), 'T', 1., ztmlres2(:,:,jn), 'T', 1. )
+ CALL lbc_lnk( 'trdmxl_trc', ztmltot2(:,:,jn), 'T', 1._wp, ztmlres2(:,:,jn), 'T', 1._wp )
DO jl = 1, jpltrd_trc
- CALL lbc_lnk( 'trdmxl_trc', ztmltrd2(:,:,jl,jn), 'T', 1. ) ! will be output in the NetCDF trends file
+ CALL lbc_lnk( 'trdmxl_trc', ztmltrd2(:,:,jl,jn), 'T', 1._wp ) ! will be output in the NetCDF trends file
END DO
ENDIF
diff --git a/tests/ADIAB_WAVE/MY_SRC/usrdef_zgr.F90 b/tests/ADIAB_WAVE/MY_SRC/usrdef_zgr.F90
index ed83f79c1c1020b8f309e0a2631672caf02665d1..892cf155ff02ce4d6f51725790f6290ec974bd1b 100644
--- a/tests/ADIAB_WAVE/MY_SRC/usrdef_zgr.F90
+++ b/tests/ADIAB_WAVE/MY_SRC/usrdef_zgr.F90
@@ -102,7 +102,7 @@ CONTAINS
DO_2D(0,0,0,0)
zhu(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji+1,jj) )
END_2D
- CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1. ) ! boundary condition: this mask the surrouding grid-points
+ CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1._wp ) ! boundary condition: this mask the surrouding grid-points
! ! ==>>> set by hand non-zero value on first/last columns & rows
DO ji = mi0(1,nn_hls), mi1(1,nn_hls) ! first row of global domain only
zhu(ji,2) = zht(ji,2)
@@ -121,7 +121,7 @@ CONTAINS
! no ocean cavities : top ocean level is ONE, except over land
! the ocean basin surrounded by land (1+nn_hls grid-points) set through lbc_lnk call
z2d(:,:) = 1._wp ! surface ocean is the 1st level
- CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! closed basin, see userdef_nam.F90
+ CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp ) ! closed basin, see userdef_nam.F90
k_top(:,:) = NINT( z2d(:,:) )
!
!
diff --git a/tests/BENCH/MY_SRC/usrdef_istate.F90 b/tests/BENCH/MY_SRC/usrdef_istate.F90
index bb35d7849422e2cc9f34673f677d4a8d5cc22af5..226b223969ecc63012c63bfc85c6c4affbf29592 100644
--- a/tests/BENCH/MY_SRC/usrdef_istate.F90
+++ b/tests/BENCH/MY_SRC/usrdef_istate.F90
@@ -82,9 +82,9 @@ CONTAINS
pu( :,:,jpk ) = 0._wp
pv( :,:,jpk ) = 0._wp
!
- CALL lbc_lnk('usrdef_istate', pts, 'T', 1. ) ! apply boundary conditions
- CALL lbc_lnk('usrdef_istate', pu, 'U', -1. ) ! apply boundary conditions
- CALL lbc_lnk('usrdef_istate', pv, 'V', -1. ) ! apply boundary conditions
+ CALL lbc_lnk('usrdef_istate', pts, 'T', 1._wp ) ! apply boundary conditions
+ CALL lbc_lnk('usrdef_istate', pu, 'U', -1._wp ) ! apply boundary conditions
+ CALL lbc_lnk('usrdef_istate', pv, 'V', -1._wp ) ! apply boundary conditions
END SUBROUTINE usr_def_istate
@@ -111,7 +111,7 @@ CONTAINS
pssh(ji,jj) = 0.1 * ( 0.5 - REAL( mig(ji,0) + (mjg(jj,0)-1) * Ni0glo, wp ) / REAL( Ni0glo * Nj0glo, wp ) )
END_2D
!
- CALL lbc_lnk('usrdef_istate', pssh, 'T', 1. ) ! apply boundary conditions
+ CALL lbc_lnk('usrdef_istate', pssh, 'T', 1._wp ) ! apply boundary conditions
!
END SUBROUTINE usr_def_istate_ssh
diff --git a/tests/BENCH/MY_SRC/usrdef_sbc.F90 b/tests/BENCH/MY_SRC/usrdef_sbc.F90
index 899c20bb3d6cfbdc8b273b08a7efb771a9fd4fc2..50ebfec3b58248a2aadcf942b8e99eb737da68f9 100644
--- a/tests/BENCH/MY_SRC/usrdef_sbc.F90
+++ b/tests/BENCH/MY_SRC/usrdef_sbc.F90
@@ -116,7 +116,7 @@ CONTAINS
vtau_ice(mi0( nn_hls+1,nn_hls):mi1( nn_hls+1,nn_hls),mj0(jpjglo-nn_hls,nn_hls):mj1(jpjglo-nn_hls,nn_hls)) = 0._wp
ENDIF
- CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'T', -1., vtau_ice, 'T', -1., ldfull = .TRUE. )
+ CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'T', -1._wp, vtau_ice, 'T', -1._wp, ldfull = .TRUE. )
#endif
!
END SUBROUTINE usrdef_sbc_ice_tau
diff --git a/tests/CANAL/MY_SRC/usrdef_istate.F90 b/tests/CANAL/MY_SRC/usrdef_istate.F90
index fca0099da7637f80e53b9425d0112546370a13cb..498f8c49962883e5d31ead56f74cbc0542642bf4 100644
--- a/tests/CANAL/MY_SRC/usrdef_istate.F90
+++ b/tests/CANAL/MY_SRC/usrdef_istate.F90
@@ -237,8 +237,8 @@ CONTAINS
!
END SELECT
!
- CALL lbc_lnk( 'usrdef_istate', pts , 'T', 1. )
- CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. )
+ CALL lbc_lnk( 'usrdef_istate', pts , 'T', 1._wp )
+ CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1._wp, pv, 'V', -1._wp )
END SUBROUTINE usr_def_istate
@@ -322,7 +322,7 @@ CONTAINS
DO_2D( 0, 0, 0, 0 )
pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * rn_uzonal * EXP( - 0.5 * gphit(ji,jj)**2 / rn_lambda**2 ) * e2t(ji,jj)
END_2D
- CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T', 1. )
+ CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T', 1._wp )
END DO
!
CASE(4) !== geostrophic zonal pulse !!st need to implement a way to separate ssh properly ==!
@@ -374,7 +374,7 @@ CONTAINS
CALL RANDOM_NUMBER(zrandom)
pssh(:,:) = pssh(:,:) + ( 0.1 * zrandom(:,:) - 0.05 )
ENDIF
- CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T', 1. )
+ CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T', 1._wp )
!
END SUBROUTINE usr_def_istate_ssh
diff --git a/tests/CANAL/MY_SRC/usrdef_zgr.F90 b/tests/CANAL/MY_SRC/usrdef_zgr.F90
index e6438852a9db61a58dc973a7a369fe7dd4663979..82d03f4329b7dc9318beab770fd1f768f8312fc7 100644
--- a/tests/CANAL/MY_SRC/usrdef_zgr.F90
+++ b/tests/CANAL/MY_SRC/usrdef_zgr.F90
@@ -199,7 +199,7 @@ CONTAINS
z2d(:,:) = REAL(jpkm1 - NINT( 0.5 * REAL(jpkm1,wp) * z2d(:,:) ), wp)
END SELECT
!
- CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! set surrounding land to zero (closed boundaries)
+ CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp ) ! set surrounding land to zero (closed boundaries)
!
k_bot(:,:) = NINT( z2d(:,:) ) ! =jpkm1 over the ocean point, =0 elsewhere
!
diff --git a/tests/DOME/MY_SRC/usrdef_zgr.F90 b/tests/DOME/MY_SRC/usrdef_zgr.F90
index 0fb7ecd0cfb333f9ec198dfd3e0f833f8bab3350..40bb48bd397c2ff01ad0ed15c626869ed284922e 100644
--- a/tests/DOME/MY_SRC/usrdef_zgr.F90
+++ b/tests/DOME/MY_SRC/usrdef_zgr.F90
@@ -103,7 +103,7 @@ CONTAINS
WHERE (gphit(:,:)>0._wp) z2d(:,:) = 0._wp
! Dig inlet:
WHERE ((gphit(:,:)>0._wp).AND.(glamt(:,:)>-50._wp).AND.(glamt(:,:)<50._wp)) z2d(:,:) = 1._wp
- CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! closed basin, see userdef_nam.F90
+ CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp ) ! closed basin, see userdef_nam.F90
k_top(:,:) = NINT( z2d(:,:) )
!
!
diff --git a/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90 b/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
index 0b4f020abb72d3849d89ea5be6bac204539047c2..0a36b514e7c25022710760c3f572110821196eb6 100644
--- a/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
+++ b/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
@@ -140,8 +140,8 @@ CONTAINS
utau_ice(ji,jj) = zrhoa * Cd_atm * wndm_ice(ji,jj) * zwndi_t
vtau_ice(ji,jj) = zrhoa * Cd_atm * wndm_ice(ji,jj) * zwndj_t
END_2D
- CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'T', -1., vtau_ice, 'T', -1., ldfull = .TRUE. )
-
+
+ CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'T', -1._wp, vtau_ice, 'T', -1._wp, ldfull = .TRUE. )
!
END SUBROUTINE usrdef_sbc_ice_tau
diff --git a/tests/LOCK_EXCHANGE/MY_SRC/usrdef_zgr.F90 b/tests/LOCK_EXCHANGE/MY_SRC/usrdef_zgr.F90
index 2c5879eba4bfe0e8a29626efe5258d9171da7b2b..5592ac6bfbca6ca0dbbd77fd087d3fe12c21350e 100644
--- a/tests/LOCK_EXCHANGE/MY_SRC/usrdef_zgr.F90
+++ b/tests/LOCK_EXCHANGE/MY_SRC/usrdef_zgr.F90
@@ -86,7 +86,7 @@ CONTAINS
! no ocean cavities : top ocean level is ONE, except over land
! the ocean basin surrounded by land (1+nn_hls grid-points) set through lbc_lnk call
z2d(:,:) = 1._wp ! surface ocean is the 1st level
- CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! closed basin, see userdef_nam.F90
+ CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp ) ! closed basin, see userdef_nam.F90
k_top(:,:) = NINT( z2d(:,:) )
!
!
diff --git a/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90 b/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90
index 9a88ea58e0a6678fe1e608294ecf1be735e7a47b..af15e507485a2e38e5367ca4f77ce6c05a4894e4 100644
--- a/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90
+++ b/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90
@@ -94,7 +94,7 @@ CONTAINS
DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls )
zhu(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji+1,jj) )
END_2D
- CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1. ) ! boundary condition: this mask the surrouding grid-points
+ CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1._wp ) ! boundary condition: this mask the surrouding grid-points
! ! ==>>> set by hand non-zero value on first/last columns & rows
!
CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d ) ! Reference z-coordinate system
@@ -105,7 +105,7 @@ CONTAINS
! no ocean cavities : top ocean level is ONE, except over land
! the ocean basin surrounded by land (1+nn_hls grid-points) set through lbc_lnk call
z2d(:,:) = 1._wp ! surface ocean is the 1st level
- CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! closed basin, see userdef_nam.F90
+ CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp ) ! closed basin, see userdef_nam.F90
k_top(:,:) = NINT( z2d(:,:) )
!
!
diff --git a/tests/SWG/MY_SRC/usrdef_zgr.F90 b/tests/SWG/MY_SRC/usrdef_zgr.F90
index 7da49f1f646f73fceb72a5d6b8e67917d523a537..9bfe7b8e061e1232e187ef2aa3c90aba160c1f72 100644
--- a/tests/SWG/MY_SRC/usrdef_zgr.F90
+++ b/tests/SWG/MY_SRC/usrdef_zgr.F90
@@ -191,7 +191,7 @@ CONTAINS
!
z2d(:,:) = REAL( jpkm1 , wp ) ! flat bottom
!
- CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! set surrounding land to zero (closed boundaries)
+ CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp ) ! set surrounding land to zero (closed boundaries)
!
zylim0 = 10000._wp ! +10km
zylim1 = 2010000._wp ! 2010km
diff --git a/tests/TSUNAMI/MY_SRC/usrdef_istate.F90 b/tests/TSUNAMI/MY_SRC/usrdef_istate.F90
index 0ca6c3492fd66504612e2c8f224d4c330d8eee84..69202f700203415cb384af1d7f562d3486fc8c27 100644
--- a/tests/TSUNAMI/MY_SRC/usrdef_istate.F90
+++ b/tests/TSUNAMI/MY_SRC/usrdef_istate.F90
@@ -103,7 +103,7 @@ CONTAINS
END_2D
!
- CALL lbc_lnk('usrdef_istate', pssh, 'T', 1. ) ! apply boundary conditions
+ CALL lbc_lnk('usrdef_istate', pssh, 'T', 1._wp ) ! apply boundary conditions
!
END SUBROUTINE usr_def_istate_ssh
diff --git a/tests/VORTEX/MY_SRC/usrdef_istate.F90 b/tests/VORTEX/MY_SRC/usrdef_istate.F90
index 6de0b19191cffc5940d9cc6763bd90c93f90ec48..7694c43c9c5ec788f0a943a2615b431b2e0a1b4c 100644
--- a/tests/VORTEX/MY_SRC/usrdef_istate.F90
+++ b/tests/VORTEX/MY_SRC/usrdef_istate.F90
@@ -122,7 +122,7 @@ CONTAINS
END DO
END_2D
!
- CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. )
+ CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1._wp, pv, 'V', -1._wp )
!
END SUBROUTINE usr_def_istate
diff --git a/tests/VORTEX/MY_SRC/usrdef_zgr.F90 b/tests/VORTEX/MY_SRC/usrdef_zgr.F90
index 06c739ed22a38ec08ce54d3f9f7ed80b18c289a0..b733421b153a4b765377650a8abbb257547fa819 100644
--- a/tests/VORTEX/MY_SRC/usrdef_zgr.F90
+++ b/tests/VORTEX/MY_SRC/usrdef_zgr.F90
@@ -187,7 +187,7 @@ CONTAINS
!
z2d(:,:) = REAL( jpkm1 , wp ) ! flat bottom
!
- CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! set surrounding land to zero (closed boundaries)
+ CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp ) ! set surrounding land to zero (closed boundaries)
!
k_bot(:,:) = NINT( z2d(:,:) ) ! =jpkm1 over the ocean point, =0 elsewhere
!
diff --git a/tests/WAD/MY_SRC/usrdef_zgr.F90 b/tests/WAD/MY_SRC/usrdef_zgr.F90
index df9970a0b289ce811f436720e5ef95cf711b0a36..3af894e85f6f7f40a00cabb3f43cfccfd95028b2 100644
--- a/tests/WAD/MY_SRC/usrdef_zgr.F90
+++ b/tests/WAD/MY_SRC/usrdef_zgr.F90
@@ -232,7 +232,7 @@ CONTAINS
DO_2D( 0, 0, 0, 0 )
zhu(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji+1,jj) )
END_2D
- CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1. ) ! boundary condition: this mask the surrounding grid-points
+ CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1._wp ) ! boundary condition: this mask the surrounding grid-points
! ! ==>>> set by hand non-zero value on first/last columns & rows
DO ji = mi0( 1,nn_hls), mi1( 1,nn_hls) ! first row of global domain only
zhu(ji,:) = zht(1,:)
@@ -245,7 +245,7 @@ CONTAINS
DO_2D( 0, 0, 0, 0 )
zhv(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji,jj+1) )
END_2D
- CALL lbc_lnk( 'usrdef_zgr', zhv, 'V', 1. ) ! boundary condition: this mask the surrounding grid-points
+ CALL lbc_lnk( 'usrdef_zgr', zhv, 'V', 1._wp ) ! boundary condition: this mask the surrounding grid-points
DO jj = mj0( 1,nn_hls), mj1( 1,nn_hls) ! first row of global domain only
zhv(:,jj) = zht(:,jj)
END DO
@@ -266,7 +266,7 @@ CONTAINS
z2d(:,mj0( 1,nn_hls):mj1( 1,nn_hls)) = 0._wp
z2d(:,mj0(jpjglo,nn_hls):mj1(jpjglo,nn_hls)) = 0._wp
- CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! closed basin, see userdef_nam.F90
+ CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp ) ! closed basin, see userdef_nam.F90
k_top(:,:) = NINT( z2d(:,:) )
!
!
@@ -301,15 +301,15 @@ CONTAINS
pe3vw(ji,jj,jk) = zwet * z1_jpkm1
END DO
END_2D
- CALL lbc_lnk( 'usrdef_zgr', pdept, 'T', 1. )
- CALL lbc_lnk( 'usrdef_zgr', pdepw, 'T', 1. )
- CALL lbc_lnk( 'usrdef_zgr', pe3t , 'T', 1. )
- CALL lbc_lnk( 'usrdef_zgr', pe3w , 'T', 1. )
- CALL lbc_lnk( 'usrdef_zgr', pe3u , 'U', 1. )
- CALL lbc_lnk( 'usrdef_zgr', pe3uw, 'U', 1. )
- CALL lbc_lnk( 'usrdef_zgr', pe3f , 'F', 1. )
- CALL lbc_lnk( 'usrdef_zgr', pe3v , 'V', 1. )
- CALL lbc_lnk( 'usrdef_zgr', pe3vw, 'V', 1. )
+ CALL lbc_lnk( 'usrdef_zgr', pdept, 'T', 1._wp )
+ CALL lbc_lnk( 'usrdef_zgr', pdepw, 'T', 1._wp )
+ CALL lbc_lnk( 'usrdef_zgr', pe3t , 'T', 1._wp )
+ CALL lbc_lnk( 'usrdef_zgr', pe3w , 'T', 1._wp )
+ CALL lbc_lnk( 'usrdef_zgr', pe3u , 'U', 1._wp )
+ CALL lbc_lnk( 'usrdef_zgr', pe3uw, 'U', 1._wp )
+ CALL lbc_lnk( 'usrdef_zgr', pe3f , 'F', 1._wp )
+ CALL lbc_lnk( 'usrdef_zgr', pe3v , 'V', 1._wp )
+ CALL lbc_lnk( 'usrdef_zgr', pe3vw, 'V', 1._wp )
WHERE( pe3t (:,:,:) == 0._wp ) pe3t (:,:,:) = 1._wp
WHERE( pe3u (:,:,:) == 0._wp ) pe3u (:,:,:) = 1._wp
WHERE( pe3v (:,:,:) == 0._wp ) pe3v (:,:,:) = 1._wp
diff --git a/tools/MISCELLANEOUS/add_wp_in_call_lbc_lnk.sh b/tools/MISCELLANEOUS/add_wp_in_call_lbc_lnk.sh
new file mode 100755
index 0000000000000000000000000000000000000000..5739b9283aa5b0ef6db462c6358f4c3f0ddc0ff5
--- /dev/null
+++ b/tools/MISCELLANEOUS/add_wp_in_call_lbc_lnk.sh
@@ -0,0 +1,15 @@
+#!/bin/bash
+#
+# this sed makes sure that we alwas specify the _wp in the calls to lbc_lnk
+#
+for ff in $( find src -type f ) */*/MY_SRC/* # for all nemo files
+do
+ # look for something like ", 'X', 1. )" and replace it by ", 'X', 1.wp )"
+ # - X can be TUVFW
+ # - the "1" can be followed by a "." or not. It must be followed by "," or ")" or "&"
+ # - any number of blanks is acceptable.
+ sed -e "s/\(, *'[TUVFW]' *, *-*1\)\.*\( *[,)&]\)/\1._wp\2/g" $ff > tmp
+ mv tmp $ff
+done
+
+