From 4e041290c37081c2c735b95310888be3b4d9825e Mon Sep 17 00:00:00 2001
From: clem <clement.rousset@locean.ipsl.fr>
Date: Tue, 4 Oct 2022 20:52:52 +0200
Subject: [PATCH] first steps of improvements of memory footprint in ice
routines
---
src/ICE/ice1d.F90 | 6 +
src/ICE/icecor.F90 | 7 +-
src/ICE/icedyn_rdgrft.F90 | 200 ++++++++-------------
src/ICE/iceitd.F90 | 118 ++++--------
src/ICE/icetab.F90 | 83 ++++++++-
src/ICE/icethd_do.F90 | 28 +--
src/ICE/icethd_zdf_bl99.F90 | 347 +++++++++++++++---------------------
7 files changed, 350 insertions(+), 439 deletions(-)
diff --git a/src/ICE/ice1d.F90 b/src/ICE/ice1d.F90
index 8c610bf49..89ac7118f 100644
--- a/src/ICE/ice1d.F90
+++ b/src/ICE/ice1d.F90
@@ -166,6 +166,9 @@ MODULE ice1D
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: a_ib_2d
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_ib_2d
+
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_i_2d
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_s_2d
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -235,6 +238,9 @@ CONTAINS
& v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) , &
& a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , v_il_2d(jpij,jpl) , &
& STAT=ierr(ii) )
+ !
+ ii = ii + 1
+ ALLOCATE( e_i_2d(jpij,nlay_i,jpl) , e_s_2d(jpij,nlay_s,jpl) , STAT=ierr(ii) )
ice1D_alloc = MAXVAL( ierr(:) )
IF( ice1D_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice1D_alloc: failed to allocate arrays.' )
diff --git a/src/ICE/icecor.F90 b/src/ICE/icecor.F90
index 0b7b53e6a..3c45c3aef 100644
--- a/src/ICE/icecor.F90
+++ b/src/ICE/icecor.F90
@@ -71,7 +71,10 @@ CONTAINS
WHERE( a_i(:,:,:) >= epsi20 ) ; h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:)
ELSEWHERE ; h_i(:,:,:) = 0._wp
END WHERE
- WHERE( h_i(:,:,:) < rn_himin ) a_i(:,:,:) = a_i(:,:,:) * h_i(:,:,:) / rn_himin
+ IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
+ WHERE( h_i(:,:,:) < rn_himin ) a_ip(:,:,:) = a_ip(:,:,:) * h_i(:,:,:) / rn_himin
+ ENDIF
+ WHERE( h_i(:,:,:) < rn_himin ) a_i(:,:,:) = a_i (:,:,:) * h_i(:,:,:) / rn_himin
!
! !-----------------------------------------------------
! ! ice concentration should not exceed amax !
@@ -83,7 +86,7 @@ CONTAINS
! !-----------------------------------------------------
! ! Rebin categories with thickness out of bounds !
! !-----------------------------------------------------
- IF ( jpl > 1 ) CALL ice_itd_reb( kt )
+ IF( jpl > 1 ) CALL ice_itd_reb( kt )
!
! !-----------------------------------------------------
IF ( nn_icesal == 2 ) THEN ! salinity must stay in bounds [Simin,Simax] !
diff --git a/src/ICE/icedyn_rdgrft.F90 b/src/ICE/icedyn_rdgrft.F90
index 0f6afe0fb..e083da9c5 100644
--- a/src/ICE/icedyn_rdgrft.F90
+++ b/src/ICE/icedyn_rdgrft.F90
@@ -56,8 +56,6 @@ MODULE icedyn_rdgrft
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hi_hrdg ! thickness of ridging ice / mean ridge thickness
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: aridge ! participating ice ridging
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: araft ! participating ice rafting
- REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_i_2d
- REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_s_2d
!
REAL(wp), PARAMETER :: hrdg_hi_min = 1.1_wp ! min ridge thickness multiplier: min(hrdg/hi)
REAL(wp), PARAMETER :: hi_hrft = 0.5_wp ! rafting multiplier: (hi/hraft)
@@ -103,7 +101,7 @@ CONTAINS
ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij), &
& apartf(jpij,0:jpl), hrmin (jpij,jpl), hraft(jpij,jpl) , aridge(jpij,jpl), &
& hrmax (jpij,jpl) , hrexp (jpij,jpl), hi_hrdg(jpij,jpl) , araft(jpij,jpl) , &
- & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
+ & STAT=ice_dyn_rdgrft_alloc )
CALL mpp_sum ( 'icedyn_rdgrft', ice_dyn_rdgrft_alloc )
IF( ice_dyn_rdgrft_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_alloc: failed to allocate arrays' )
@@ -573,9 +571,9 @@ CONTAINS
REAL(wp), DIMENSION(jpij) :: airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg ! area etc of new ridges
REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft ! area etc of rafted ice
!
- REAL(wp), DIMENSION(jpij) :: ersw ! enth of water trapped into ridges
- REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2
- REAL(wp), DIMENSION(jpij) :: z1_ai ! 1 / a
+ REAL(wp) :: ersw ! enthalpy of water trapped into ridges
+ REAL(wp) :: zswitch, fvol ! new ridge volume going to jl2
+ REAL(wp) :: z1_ai ! 1 / a
REAL(wp), DIMENSION(jpij) :: zvti ! sum(v_i)
!
REAL(wp), DIMENSION(jpij,nlay_s) :: esrft ! snow energy of rafting ice
@@ -612,8 +610,8 @@ CONTAINS
IF( ll_shift(ji) ) THEN ! only if ice is ridging
- IF( a_i_2d(ji,jl1) > epsi10 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
- ELSE ; z1_ai(ji) = 0._wp
+ IF( a_i_2d(ji,jl1) > epsi10 ) THEN ; z1_ai = 1._wp / a_i_2d(ji,jl1)
+ ELSE ; z1_ai = 0._wp
ENDIF
! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
@@ -624,15 +622,15 @@ CONTAINS
airft2(ji) = airft1 * hi_hrft
! ridging /rafting fractions
- afrdg = airdg1 * z1_ai(ji)
- afrft = airft1 * z1_ai(ji)
+ afrdg = airdg1 * z1_ai
+ afrft = airft1 * z1_ai
! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
IF ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg ! v <= 10m then porosity = rn_porordg
ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp ! v >= 20m then porosity = 0
ELSE ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0
ENDIF
- ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji) ! clem: if sst>0, then ersw <0 (is that possible?)
+ ersw = -rhoi * vsw * rcp * sst_1d(ji) ! clem: if sst>0, then ersw <0 (is that possible?)
! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
virdg1 = v_i_2d (ji,jl1) * afrdg
@@ -661,16 +659,29 @@ CONTAINS
vlrft (ji) = v_il_2d(ji,jl1) * afrft
ENDIF
ENDIF
+
+ DO jk = 1, nlay_s
+ esrdg(ji,jk) = e_s_2d (ji,jk,jl1) * afrdg
+ esrft(ji,jk) = e_s_2d (ji,jk,jl1) * afrft
+ END DO
+ DO jk = 1, nlay_i
+ eirdg(ji,jk) = e_i_2d (ji,jk,jl1) * afrdg + ersw * r1_nlay_i
+ eirft(ji,jk) = e_i_2d (ji,jk,jl1) * afrft
+ END DO
! Ice-ocean exchanges associated with ice porosity
wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoi * r1_Dt_ice ! increase in ice volume due to seawater frozen in voids
sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_Dt_ice
- hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_Dt_ice ! > 0 [W.m-2]
+ hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw * r1_Dt_ice ! > 0 [W.m-2]
! Put the snow lost by ridging into the ocean
! Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow.
wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhos * vsrdg(ji) * ( 1._wp - rn_fsnwrdg ) & ! fresh water source for ocean
& + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice
+ DO jk = 1, nlay_s
+ hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg ) & ! heat sink for ocean (<0, W.m-2)
+ & - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice
+ END DO
! virtual salt flux to keep salinity constant
IF( nn_icesal /= 2 ) THEN
@@ -693,49 +704,18 @@ CONTAINS
v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - vlrdg(ji) - vlrft(ji)
ENDIF
ENDIF
+ DO jk = 1, nlay_s
+ e_s_2d(ji,jk,jl1) = e_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )
+ END DO
+ DO jk = 1, nlay_i
+ e_i_2d(ji,jk,jl1) = e_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )
+ END DO
+
ENDIF
-
+
END DO ! ji
-
- ! special loop for e_s because of layers jk
- DO jk = 1, nlay_s
- DO ji = 1, npti
- IF( ll_shift(ji) ) THEN
- ! Compute ridging /rafting fractions
- afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
- afrft = araft (ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
- ! Compute ridging /rafting ice and new ridges for es
- esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
- esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
- ! Put the snow lost by ridging into the ocean
- hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg ) & ! heat sink for ocean (<0, W.m-2)
- & - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice
- !
- ! Remove energy of new ridge to each category jl1
- !-------------------------------------------------
- ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )
- ENDIF
- END DO
- END DO
-
- ! special loop for e_i because of layers jk
- DO jk = 1, nlay_i
- DO ji = 1, npti
- IF( ll_shift(ji) ) THEN
- ! Compute ridging /rafting fractions
- afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
- afrft = araft (ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
- ! Compute ridging ice and new ridges for ei
- eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
- eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
- !
- ! Remove energy of new ridge to each category jl1
- !-------------------------------------------------
- ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )
- ENDIF
- END DO
- END DO
-
+
+
! 3) compute categories in which ice is added (jl2)
!--------------------------------------------------
itest_rdg(1:npti) = 0
@@ -752,13 +732,13 @@ CONTAINS
IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
hR = MIN( hrmax(ji,jl1), hi_max(jl2) )
- farea = ( hR - hL ) / ( hrmax(ji,jl1) - hrmin(ji,jl1) )
- fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
+ farea = ( hR - hL ) / ( hrmax(ji,jl1) - hrmin(ji,jl1) )
+ fvol = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
!
itest_rdg(ji) = 1 ! test for conservation
ELSE
- farea = 0._wp
- fvol(ji) = 0._wp
+ farea = 0._wp
+ fvol = 0._wp
ENDIF
!
ELSEIF( ln_distf_exp ) THEN ! Lipscomb et al. (2007) exponential formulation
@@ -766,24 +746,24 @@ CONTAINS
IF( jl2 < jpl ) THEN
!
IF( hrmin(ji,jl1) <= hi_max(jl2) ) THEN
- hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
- hR = hi_max(jl2)
- expL = EXP( -( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
- expR = EXP( -( hR - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
- farea = expL - expR
- fvol(ji) = ( ( hL + hrexp(ji,jl1) ) * expL &
- - ( hR + hrexp(ji,jl1) ) * expR ) / ( hrmin(ji,jl1) + hrexp(ji,jl1) )
+ hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
+ hR = hi_max(jl2)
+ expL = EXP( -( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
+ expR = EXP( -( hR - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
+ farea = expL - expR
+ fvol = ( ( hL + hrexp(ji,jl1) ) * expL &
+ & - ( hR + hrexp(ji,jl1) ) * expR ) / ( hrmin(ji,jl1) + hrexp(ji,jl1) )
ELSE
- farea = 0._wp
- fvol(ji) = 0._wp
+ farea = 0._wp
+ fvol = 0._wp
END IF
!
ELSE ! jl2 = jpl
!
- hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
- expL = EXP(-( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
- farea = expL
- fvol(ji) = ( hL + hrexp(ji,jl1) ) * expL / ( hrmin(ji,jl1) + hrexp(ji,jl1) )
+ hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
+ expL = EXP(-( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) )
+ farea = expL
+ fvol = ( hL + hrexp(ji,jl1) ) * expL / ( hrmin(ji,jl1) + hrexp(ji,jl1) )
!
END IF ! jl2 < jpl
!
@@ -793,59 +773,49 @@ CONTAINS
! Compute the fraction of rafted ice area and volume going to thickness category jl2
IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) > hi_max(jl2-1) ) THEN
- zswitch(ji) = 1._wp
+ zswitch = 1._wp
!
itest_rft(ji) = 1 ! test for conservation
ELSE
- zswitch(ji) = 0._wp
+ zswitch = 0._wp
ENDIF
!
! Patch to ensure perfect conservation if ice thickness goes mad
! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas)
! Then ice volume is removed from one category but the ridging/rafting scheme
! does not know where to move it, leading to a conservation issue.
- IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN ; farea = 1._wp ; fvol(ji) = 1._wp ; ENDIF
- IF( itest_rft(ji) == 0 .AND. jl2 == jpl ) zswitch(ji) = 1._wp
+ IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN ; farea = 1._wp ; fvol = 1._wp ; ENDIF
+ IF( itest_rft(ji) == 0 .AND. jl2 == jpl ) zswitch = 1._wp
!
! Add area, volume of new ridge to category jl2
!----------------------------------------------
- a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea + airft2(ji) * zswitch(ji) )
- oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea + oirft2(ji) * zswitch(ji) )
- v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
- sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
- v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji) + &
- & vsrft (ji) * rn_fsnwrft * zswitch(ji) )
+ a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea + airft2(ji) * zswitch )
+ oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea + oirft2(ji) * zswitch )
+ v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol + virft (ji) * zswitch )
+ sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol + sirft (ji) * zswitch )
+ v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol + &
+ & vsrft (ji) * rn_fsnwrft * zswitch )
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
- v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + ( vprdg (ji) * rn_fpndrdg * fvol (ji) &
- & + vprft (ji) * rn_fpndrft * zswitch(ji) )
- a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + ( aprdg2(ji) * rn_fpndrdg * farea &
- & + aprft2(ji) * rn_fpndrft * zswitch(ji) )
+ v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + ( vprdg (ji) * rn_fpndrdg * fvol &
+ & + vprft (ji) * rn_fpndrft * zswitch )
+ a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + ( aprdg2(ji) * rn_fpndrdg * farea &
+ & + aprft2(ji) * rn_fpndrft * zswitch )
IF ( ln_pnd_lids ) THEN
- v_il_2d (ji,jl2) = v_il_2d(ji,jl2) + ( vlrdg(ji) * rn_fpndrdg * fvol (ji) &
- & + vlrft(ji) * rn_fpndrft * zswitch(ji) )
+ v_il_2d (ji,jl2) = v_il_2d(ji,jl2) + ( vlrdg(ji) * rn_fpndrdg * fvol &
+ & + vlrft(ji) * rn_fpndrft * zswitch )
ENDIF
ENDIF
-
+ DO jk = 1, nlay_s
+ e_s_2d(ji,jk,jl2) = e_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol + &
+ & esrft(ji,jk) * rn_fsnwrft * zswitch )
+ END DO
+ DO jk = 1, nlay_i
+ e_i_2d(ji,jk,jl2) = e_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol + eirft(ji,jk) * zswitch
+ END DO
+
ENDIF
END DO
- ! Add snow energy of new ridge to category jl2
- !---------------------------------------------
- DO jk = 1, nlay_s
- DO ji = 1, npti
- IF( ll_shift(ji) ) &
- & ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji) + &
- & esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
- END DO
- END DO
- ! Add ice energy of new ridge to category jl2
- !--------------------------------------------
- DO jk = 1, nlay_i
- DO ji = 1, npti
- IF( ll_shift(ji) ) &
- & ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)
- END DO
- END DO
!
END DO ! jl2
!
@@ -854,7 +824,7 @@ CONTAINS
! roundoff errors
!----------------
! In case ridging/rafting lead to very small negative values (sometimes it happens)
- CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, ze_s_2d, ze_i_2d )
+ CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, e_s_2d, e_i_2d )
!
END SUBROUTINE rdgrft_shift
@@ -1036,14 +1006,8 @@ CONTAINS
CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
- DO jl = 1, jpl
- DO jk = 1, nlay_s
- CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
- END DO
- DO jk = 1, nlay_i
- CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
- END DO
- END DO
+ CALL tab_4d_3d( npti, nptidx(1:npti), e_s_2d (1:npti,:,:) , e_s )
+ CALL tab_4d_3d( npti, nptidx(1:npti), e_i_2d (1:npti,:,:) , e_i )
CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) )
CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) )
CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) )
@@ -1063,14 +1027,8 @@ CONTAINS
CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
- DO jl = 1, jpl
- DO jk = 1, nlay_s
- CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
- END DO
- DO jk = 1, nlay_i
- CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
- END DO
- END DO
+ CALL tab_3d_4d( npti, nptidx(1:npti), e_s_2d (1:npti,:,:) , e_s )
+ CALL tab_3d_4d( npti, nptidx(1:npti), e_i_2d (1:npti,:,:) , e_i )
CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) )
CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) )
CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) )
diff --git a/src/ICE/iceitd.F90 b/src/ICE/iceitd.F90
index caa37df09..e5ed9a341 100644
--- a/src/ICE/iceitd.F90
+++ b/src/ICE/iceitd.F90
@@ -306,25 +306,6 @@ CONTAINS
! 6) Shift ice between categories
!----------------------------------------------------------------------------------------------
CALL itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )
-
- !----------------------------------------------------------------------------------------------
- ! 7) Make sure h_i >= minimum ice thickness hi_min
- !----------------------------------------------------------------------------------------------
- CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,1) )
- CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,1) )
- CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,1) )
- !
- DO ji = 1, npti
- IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
- a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin
- IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
- h_i_1d(ji) = rn_himin
- ENDIF
- END DO
- !
- CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,1) )
- CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,1) )
- CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,1) )
!
ENDIF
!
@@ -415,11 +396,9 @@ CONTAINS
!
INTEGER :: ji, jl, jk ! dummy loop indices
INTEGER :: jl2, jl1 ! local integers
- REAL(wp) :: ztrans ! ice/snow transferred
- REAL(wp), DIMENSION(jpij) :: zworka, zworkv ! workspace
+ REAL(wp) :: zworka, zworkv, ztrans ! ice/snow transferred
+ REAL(wp), DIMENSION(jpij) :: ztmp ! workspace
REAL(wp), DIMENSION(jpij,jpl) :: zaTsfn ! - -
- REAL(wp), DIMENSION(jpij,nlay_i,jpl) :: ze_i_2d
- REAL(wp), DIMENSION(jpij,nlay_s,jpl) :: ze_s_2d
!!------------------------------------------------------------------
CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i )
@@ -432,14 +411,8 @@ CONTAINS
CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il )
CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
- DO jl = 1, jpl
- DO jk = 1, nlay_s
- CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
- END DO
- DO jk = 1, nlay_i
- CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
- END DO
- END DO
+ CALL tab_4d_3d( npti, nptidx(1:npti), e_s_2d (1:npti,:,:) , e_s )
+ CALL tab_4d_3d( npti, nptidx(1:npti), e_i_2d (1:npti,:,:) , e_i )
! to correct roundoff errors on a_i
CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d )
@@ -466,11 +439,11 @@ CONTAINS
ELSE ; jl2 = jl
ENDIF
!
- IF( v_i_2d(ji,jl1) >= epsi10 ) THEN ; zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
- ELSE ; zworkv(ji) = 0._wp
+ IF( v_i_2d(ji,jl1) >= epsi10 ) THEN ; zworkv = pdvice(ji,jl) / v_i_2d(ji,jl1)
+ ELSE ; zworkv = 0._wp
ENDIF
- IF( a_i_2d(ji,jl1) >= epsi10 ) THEN ; zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
- ELSE ; zworka(ji) = 0._wp
+ IF( a_i_2d(ji,jl1) >= epsi10 ) THEN ; zworka = pdaice(ji,jl) / a_i_2d(ji,jl1)
+ ELSE ; zworka = 0._wp
ENDIF
!
a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) ! Ice areas
@@ -479,71 +452,50 @@ CONTAINS
v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl) ! Ice volumes
v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
!
- ztrans = v_s_2d(ji,jl1) * zworkv(ji) ! Snow volumes
+ ztrans = v_s_2d(ji,jl1) * zworkv ! Snow volumes
v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans
!
- ztrans = oa_i_2d(ji,jl1) * zworka(ji) ! Ice age
+ ztrans = oa_i_2d(ji,jl1) * zworka ! Ice age
oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - ztrans
oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ztrans
!
- ztrans = sv_i_2d(ji,jl1) * zworkv(ji) ! Ice salinity
+ ztrans = sv_i_2d(ji,jl1) * zworkv ! Ice salinity
sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - ztrans
sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ztrans
!
- ztrans = zaTsfn(ji,jl1) * zworka(ji) ! Surface temperature
+ ztrans = zaTsfn(ji,jl1) * zworka ! Surface temperature
zaTsfn(ji,jl1) = zaTsfn(ji,jl1) - ztrans
zaTsfn(ji,jl2) = zaTsfn(ji,jl2) + ztrans
!
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
- ztrans = a_ip_2d(ji,jl1) * zworka(ji) ! Pond fraction
+ ztrans = a_ip_2d(ji,jl1) * zworka ! Pond fraction
a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
!
- ztrans = v_ip_2d(ji,jl1) * zworkv(ji) ! Pond volume
+ ztrans = v_ip_2d(ji,jl1) * zworkv ! Pond volume
v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
!
IF ( ln_pnd_lids ) THEN ! Pond lid volume
- ztrans = v_il_2d(ji,jl1) * zworkv(ji)
+ ztrans = v_il_2d(ji,jl1) * zworkv
v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans
v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans
ENDIF
ENDIF
!
- ENDIF ! jl1 >0
- END DO
- !
- DO jk = 1, nlay_s !--- Snow heat content
- DO ji = 1, npti
- !
- jl1 = kdonor(ji,jl)
+ DO jk = 1, nlay_s ! Snow heat content
+ ztrans = e_s_2d(ji,jk,jl1) * zworkv
+ e_s_2d(ji,jk,jl1) = e_s_2d(ji,jk,jl1) - ztrans
+ e_s_2d(ji,jk,jl2) = e_s_2d(ji,jk,jl2) + ztrans
+ ENDDO
+ DO jk = 1, nlay_i ! Ice heat content
+ ztrans = e_i_2d(ji,jk,jl1) * zworkv
+ e_i_2d(ji,jk,jl1) = e_i_2d(ji,jk,jl1) - ztrans
+ e_i_2d(ji,jk,jl2) = e_i_2d(ji,jk,jl2) + ztrans
+ ENDDO
!
- IF( jl1 > 0 ) THEN
- IF(jl1 == jl) THEN ; jl2 = jl+1
- ELSE ; jl2 = jl
- ENDIF
- ztrans = ze_s_2d(ji,jk,jl1) * zworkv(ji)
- ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) - ztrans
- ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ztrans
- ENDIF
- END DO
- END DO
- !
- DO jk = 1, nlay_i !--- Ice heat content
- DO ji = 1, npti
- !
- jl1 = kdonor(ji,jl)
- !
- IF( jl1 > 0 ) THEN
- IF(jl1 == jl) THEN ; jl2 = jl+1
- ELSE ; jl2 = jl
- ENDIF
- ztrans = ze_i_2d(ji,jk,jl1) * zworkv(ji)
- ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) - ztrans
- ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + ztrans
- ENDIF
- END DO
+ ENDIF ! jl1 >0
END DO
!
END DO ! boundaries, 1 to jpl-1
@@ -553,13 +505,13 @@ CONTAINS
!-------------------
! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)
! because of truncation error ( i.e. 1. - 1. /= 0 )
- CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, ze_s_2d, ze_i_2d )
+ CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, e_s_2d, e_i_2d )
! at_i must be <= rn_amax
- zworka(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
+ ztmp(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
DO jl = 1, jpl
- WHERE( zworka(1:npti) > rn_amax_1d(1:npti) ) &
- & a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti)
+ WHERE( ztmp(1:npti) > rn_amax_1d(1:npti) ) &
+ & a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / ztmp(1:npti)
END DO
!-------------------------------------------------------------------------------
@@ -587,14 +539,8 @@ CONTAINS
CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il )
CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
- DO jl = 1, jpl
- DO jk = 1, nlay_s
- CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
- END DO
- DO jk = 1, nlay_i
- CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
- END DO
- END DO
+ CALL tab_3d_4d( npti, nptidx(1:npti), e_s_2d (1:npti,:,:) , e_s )
+ CALL tab_3d_4d( npti, nptidx(1:npti), e_i_2d (1:npti,:,:) , e_i )
!
END SUBROUTINE itd_shiftice
diff --git a/src/ICE/icetab.F90 b/src/ICE/icetab.F90
index 02d3501ea..5db21b305 100644
--- a/src/ICE/icetab.F90
+++ b/src/ICE/icetab.F90
@@ -20,8 +20,10 @@ MODULE icetab
IMPLICIT NONE
PRIVATE
+ PUBLIC tab_4d_3d
PUBLIC tab_3d_2d
PUBLIC tab_2d_1d
+ PUBLIC tab_3d_4d
PUBLIC tab_2d_3d
PUBLIC tab_1d_2d
@@ -32,6 +34,39 @@ MODULE icetab
!!----------------------------------------------------------------------
CONTAINS
+ SUBROUTINE tab_4d_3d( ndim1d, tab_ind, tab1d, tab2d )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE tab_2d_1d ***
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: ndim1d ! 1d size
+ INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
+ REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: tab2d ! input 2D field
+ REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: tab1d ! output 1D field
+ !
+ INTEGER :: ipi, ipj, ipk, ji0, jj0, jk, jl, jn, jid, jjd
+ !!----------------------------------------------------------------------
+ ipi = SIZE(tab2d,1) ! 1st dimension
+ ipj = SIZE(tab2d,2) ! 2nd dimension
+ ipk = SIZE(tab2d,3) ! 3d dimension
+ !
+ IF( ipi == jpi .AND. ipj == jpj ) THEN ! full arrays then no need to change index jid and jjd
+ ji0 = 0 ; jj0 = 0
+ ELSE ! reduced arrays then need to shift index by nn_hls
+ ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls
+ ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls)
+ !
+ DO jn = 1, ndim1d
+ jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
+ jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
+ DO jl = 1, jpl
+ DO jk = 1, ipk
+ tab1d(jn,jk,jl) = tab2d(jid,jjd,jk,jl)
+ END DO
+ END DO
+ END DO
+ !
+ END SUBROUTINE tab_4d_3d
+
SUBROUTINE tab_3d_2d( ndim1d, tab_ind, tab1d, tab2d )
!!----------------------------------------------------------------------
!! *** ROUTINE tab_2d_1d ***
@@ -52,10 +87,10 @@ CONTAINS
ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls
ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls)
!
- DO jl = 1, jpl
- DO jn = 1, ndim1d
- jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
- jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
+ DO jn = 1, ndim1d
+ jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
+ jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
+ DO jl = 1, jpl
tab1d(jn,jl) = tab2d(jid,jjd,jl)
END DO
END DO
@@ -89,6 +124,38 @@ CONTAINS
END DO
END SUBROUTINE tab_2d_1d
+ SUBROUTINE tab_3d_4d( ndim1d, tab_ind, tab1d, tab2d )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE tab_2d_1d ***
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: ndim1d ! 1D size
+ INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
+ REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: tab1d ! input 1D field
+ REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: tab2d ! output 2D field
+ !
+ INTEGER :: ipi, ipj, ipk, ji0, jj0, jk, jl, jn, jid, jjd
+ !!----------------------------------------------------------------------
+ ipi = SIZE(tab2d,1) ! 1st dimension
+ ipj = SIZE(tab2d,2) ! 2nd dimension
+ ipk = SIZE(tab2d,3) ! 3d dimension
+ !
+ IF( ipi == jpi .AND. ipj == jpj ) THEN ! full arrays then no need to change index jid and jjd
+ ji0 = 0 ; jj0 = 0
+ ELSE ! reduced arrays then need to shift index by nn_hls
+ ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls
+ ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls)
+ !
+ DO jn = 1, ndim1d
+ jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
+ jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
+ DO jl = 1, jpl
+ DO jk = 1, ipk
+ tab2d(jid,jjd,jk,jl) = tab1d(jn,jk,jl)
+ END DO
+ END DO
+ END DO
+ !
+ END SUBROUTINE tab_3d_4d
SUBROUTINE tab_2d_3d( ndim1d, tab_ind, tab1d, tab2d )
!!----------------------------------------------------------------------
@@ -110,10 +177,10 @@ CONTAINS
ji0 = nn_hls ; jj0 = nn_hls ! since tab2d is shifted by nn_hls
ENDIF ! (i.e. from hls+1:jpi-hls to 1:jpi-2*hls)
!
- DO jl = 1, jpl
- DO jn = 1, ndim1d
- jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
- jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
+ DO jn = 1, ndim1d
+ jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 - ji0
+ jjd = ( tab_ind(jn) - 1 ) / jpi + 1 - jj0
+ DO jl = 1, jpl
tab2d(jid,jjd,jl) = tab1d(jn,jl)
END DO
END DO
diff --git a/src/ICE/icethd_do.F90 b/src/ICE/icethd_do.F90
index 690af0e21..8af4064e8 100644
--- a/src/ICE/icethd_do.F90
+++ b/src/ICE/icethd_do.F90
@@ -97,8 +97,6 @@ CONTAINS
REAL(wp), DIMENSION(jpij,jpl) :: zv_b ! old volume of ice in category jl
REAL(wp), DIMENSION(jpij,jpl) :: za_b ! old area of ice in category jl
!
- REAL(wp), DIMENSION(jpij,nlay_i,jpl) :: ze_i_2d !: 1-D version of e_i
- !
!!-----------------------------------------------------------------------!
IF( ln_icediachk ) CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
@@ -126,11 +124,7 @@ CONTAINS
CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
- DO jl = 1, jpl
- DO jk = 1, nlay_i
- CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
- END DO
- END DO
+ CALL tab_4d_3d( npti, nptidx(1:npti), e_i_2d (1:npti,:,:) , e_i )
CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d (1:npti) , qlead )
CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d (1:npti) , t_bo )
CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d (1:npti) , sfx_opw )
@@ -147,9 +141,9 @@ CONTAINS
DO jl = 1, jpl
DO jk = 1, nlay_i
WHERE( v_i_2d(1:npti,jl) > 0._wp )
- ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) / v_i_2d(1:npti,jl) * REAL( nlay_i )
+ e_i_2d(1:npti,jk,jl) = e_i_2d(1:npti,jk,jl) / v_i_2d(1:npti,jl) * REAL( nlay_i )
ELSEWHERE
- ze_i_2d(1:npti,jk,jl) = 0._wp
+ e_i_2d(1:npti,jk,jl) = 0._wp
END WHERE
END DO
END DO
@@ -261,8 +255,8 @@ CONTAINS
DO ji = 1, npti
jl = jcat(ji)
rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
- ze_i_2d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + &
- & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_2d(ji,jk,jl) * zv_b(ji,jl) ) &
+ e_i_2d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + &
+ & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + e_i_2d(ji,jk,jl) * zv_b(ji,jl) ) &
& * rswitch / MAX( v_i_2d(ji,jl), epsi20 )
END DO
END DO
@@ -276,7 +270,7 @@ CONTAINS
DO jk = 1, nlay_i
DO ji = 1, npti
h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i
- eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk)
+ eh_i_old(ji,jk) = e_i_2d(ji,jk,jl) * h_i_old(ji,jk)
END DO
END DO
@@ -291,7 +285,7 @@ CONTAINS
eh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
END DO
! --- Ice enthalpy remapping --- !
- CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) )
+ CALL ice_thd_ent( e_i_2d(1:npti,:,jl) )
END DO
! --- Update salinity --- !
@@ -304,7 +298,7 @@ CONTAINS
! Change units for e_i
DO jl = 1, jpl
DO jk = 1, nlay_i
- ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) * v_i_2d(1:npti,jl) * r1_nlay_i
+ e_i_2d(1:npti,jk,jl) = e_i_2d(1:npti,jk,jl) * v_i_2d(1:npti,jl) * r1_nlay_i
END DO
END DO
@@ -312,11 +306,7 @@ CONTAINS
CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
- DO jl = 1, jpl
- DO jk = 1, nlay_i
- CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
- END DO
- END DO
+ CALL tab_3d_4d( npti, nptidx(1:npti), e_i_2d (1:npti,:,:) , e_i )
CALL tab_1d_2d( npti, nptidx(1:npti), sfx_opw_1d(1:npti), sfx_opw )
CALL tab_1d_2d( npti, nptidx(1:npti), wfx_opw_1d(1:npti), wfx_opw )
CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d(1:npti), hfx_thd )
diff --git a/src/ICE/icethd_zdf_bl99.F90 b/src/ICE/icethd_zdf_bl99.F90
index 493435f3f..9bae2473e 100644
--- a/src/ICE/icethd_zdf_bl99.F90
+++ b/src/ICE/icethd_zdf_bl99.F90
@@ -105,22 +105,21 @@ CONTAINS
REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness
REAL(wp), DIMENSION(jpij) :: zh_s, z1_h_s ! snow layer thickness
REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface
- REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function
REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function
- !
- REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice
+ REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow
+ REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow
+ REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice
+ REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice
+ REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity
+ REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i_cp ! copy
REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice
REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow
+ !
+ REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function
REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Temporary temperature in the ice to check the convergence
REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence
- REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity
- REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i_cp ! copy
- REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice
- REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice
REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice
REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice
- REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow
- REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow
REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow
REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow
REAL(wp), DIMENSION(jpij) :: zkappa_comb ! Combined snow and ice surface conductivity
@@ -195,7 +194,6 @@ CONTAINS
! Store initial temperatures and non solar heat fluxes
IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
ztsub (1:npti) = t_su_1d(1:npti) ! surface temperature at iteration n-1
- ztsuold (1:npti) = t_su_1d(1:npti) ! surface temperature initial value
t_su_1d (1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! required to leave the choice between melting or not
zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux
zqns_ice_b (1:npti) = qns_ice_1d(1:npti) ! store previous qns_ice_1d value
@@ -209,19 +207,18 @@ CONTAINS
! 2) Radiation
!-------------
! --- Transmission/absorption of solar radiation in the ice --- !
- zradtr_s(1:npti,0) = qtr_ice_top_1d(1:npti)
- DO jk = 1, nlay_s
- DO ji = 1, npti
+ DO ji = 1, npti
+ !
+ zradtr_s(ji,0) = qtr_ice_top_1d(ji)
+ DO jk = 1, nlay_s
! ! radiation transmitted below the layer-th snow layer
zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s(ji) * MAX( 0._wp, zh_s(ji) * REAL(jk) - zhs_ssl ) )
! ! radiation absorbed by the layer-th snow layer
zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
END DO
- END DO
- !
- zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * za_s_fra(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - za_s_fra(1:npti) )
- DO jk = 1, nlay_i
- DO ji = 1, npti
+ !
+ zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * za_s_fra(ji) + qtr_ice_top_1d(ji) * ( 1._wp - za_s_fra(ji) )
+ DO jk = 1, nlay_i
! ! radiation transmitted below the layer-th ice layer
zradtr_i(ji,jk) = za_s_fra(ji) * zradtr_s(ji,nlay_s) & ! part covered by snow
& * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zh_min ) ) &
@@ -230,9 +227,11 @@ CONTAINS
! ! radiation absorbed by the layer-th ice layer
zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
END DO
+ !
+ qtr_ice_bot_1d(ji) = zradtr_i(ji,nlay_i) ! record radiation transmitted below the ice
+ !
END DO
!
- qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice
!
iconv = 0 ! number of iterations
!
@@ -257,9 +256,7 @@ CONTAINS
DO ji = 1, npti
ztcond_i_cp(ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
ztcond_i_cp(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )
- END DO
- DO jk = 1, nlay_i-1
- DO ji = 1, npti
+ DO jk = 1, nlay_i-1
ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &
& MIN( -epsi10, 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 )
END DO
@@ -272,9 +269,7 @@ CONTAINS
& - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) &
& - 0.011_wp * ( t_bo_1d(ji) - rt0 )
- END DO
- DO jk = 1, nlay_i-1
- DO ji = 1, npti
+ DO jk = 1, nlay_i-1
ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &
& MIN( -epsi10, 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 ) &
& - 0.011_wp * ( 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 )
@@ -308,119 +303,98 @@ CONTAINS
!
ENDIF
!
- !-----------------
- ! 4) kappa factors
- !-----------------
- !--- Snow
- ! Variable used after iterations
- ! Value must be frozen after convergence for MPP independance reason
- DO jk = 0, nlay_s-1
- DO ji = 1, npti
- IF ( .NOT. l_T_converged(ji) ) &
+ DO ji = 1, npti
+ !-----------------
+ ! 4) kappa factors
+ !-----------------
+ IF ( .NOT. l_T_converged(ji) ) THEN
+ !--- Snow
+ ! Variable used after iterations
+ ! Value must be frozen after convergence for MPP independance reason
+ DO jk = 0, nlay_s-1
zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
- END DO
- END DO
- DO ji = 1, npti ! Snow-ice interface
- IF ( .NOT. l_T_converged(ji) ) &
- zkappa_s(ji,nlay_s) = isnow(ji) * zghe(ji) * rn_cnd_s * ztcond_i(ji,0) &
+ END DO
+ zkappa_s(ji,nlay_s) = isnow(ji) * zghe(ji) * rn_cnd_s * ztcond_i(ji,0) & ! Snow-ice interface
& / ( 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) )
- END DO
-
- !--- Ice
- ! Variable used after iterations
- ! Value must be frozen after convergence for MPP independance reason
- DO jk = 0, nlay_i
- DO ji = 1, npti
- IF ( .NOT. l_T_converged(ji) ) &
+ !
+ !--- Ice
+ ! Variable used after iterations
+ ! Value must be frozen after convergence for MPP independance reason
+ DO jk = 0, nlay_i
zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
- END DO
- END DO
- DO ji = 1, npti ! Snow-ice interface
- IF ( .NOT. l_T_converged(ji) ) THEN
+ END DO
! Calculate combined surface snow and ice conductivity to pass through the coupler (met-office)
zkappa_comb(ji) = isnow_comb(ji) * zkappa_s(ji,0) + ( 1._wp - isnow_comb(ji) ) * zkappa_i(ji,0)
! If there is snow then use the same snow-ice interface conductivity for the top layer of ice
- IF( h_s_1d(ji) > 0._wp ) zkappa_i(ji,0) = zkappa_s(ji,nlay_s)
- ENDIF
- END DO
- !
- !--------------------------------------
- ! 5) Sea ice specific heat, eta factors
- !--------------------------------------
- DO jk = 1, nlay_i
- DO ji = 1, npti
+ IF( h_s_1d(ji) > 0._wp ) zkappa_i(ji,0) = zkappa_s(ji,nlay_s) ! Snow-ice interface
+ ENDIF
+ !
+ !--------------------------------------
+ ! 5) Sea ice specific heat, eta factors
+ !--------------------------------------
+ DO jk = 1, nlay_i
zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
zeta_i(ji,jk) = rDt_ice * r1_rhoi * z1_h_i(ji) / zcpi
END DO
- END DO
-
- DO jk = 1, nlay_s
- DO ji = 1, npti
+ DO jk = 1, nlay_s
zeta_s(ji,jk) = rDt_ice * r1_rhos * r1_rcpi * z1_h_s(ji)
END DO
END DO
!
- !----------------------------------------!
- ! !
- ! Conduction flux is off or emulated !
- ! !
- !----------------------------------------!
!
IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
+ !----------------------------------------!
+ ! !
+ ! Conduction flux is off or emulated !
+ ! !
+ !----------------------------------------!
!
! ==> The original BL99 temperature computation is used
! (with qsr_ice, qns_ice and dqns_ice as inputs)
!
- !----------------------------
- ! 6) surface flux computation
- !----------------------------
- ! update of the non solar flux according to the update in T_su
DO ji = 1, npti
+ !----------------------------
+ ! 6) surface flux computation
+ !----------------------------
+ ! update of the non solar flux according to the update in T_su
! Variable used after iterations
! Value must be frozen after convergence for MPP independance reason
IF ( .NOT. l_T_converged(ji) ) &
qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
- END DO
-
- DO ji = 1, npti
+ !
zfnet(ji) = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar
- END DO
- !
- !----------------------------
- ! 7) tridiagonal system terms
- !----------------------------
- ! layer denotes the number of the layer in the snow or in the ice
- ! jm denotes the reference number of the equation in the tridiagonal
- ! system, terms of tridiagonal system are indexed as following :
- ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
-
- ! ice interior terms (top equation has the same form as the others)
- ztrid (1:npti,:,:) = 0._wp
- zindterm(1:npti,:) = 0._wp
- zindtbis(1:npti,:) = 0._wp
- zdiagbis(1:npti,:) = 0._wp
-
- DO jm = nlay_s + 2, nlay_s + nlay_i
- DO ji = 1, npti
+ !
+ !
+ !----------------------------
+ ! 7) tridiagonal system terms
+ !----------------------------
+ ! layer denotes the number of the layer in the snow or in the ice
+ ! jm denotes the reference number of the equation in the tridiagonal
+ ! system, terms of tridiagonal system are indexed as following :
+ ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
+
+ ! ice interior terms (top equation has the same form as the others)
+ ztrid (ji,:,:) = 0._wp
+ zindterm(ji,:) = 0._wp
+ zindtbis(ji,:) = 0._wp
+ zdiagbis(ji,:) = 0._wp
+
+ DO jm = nlay_s + 2, nlay_s + nlay_i
jk = jm - nlay_s - 1
ztrid (ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
ztrid (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
ztrid (ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk)
zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
END DO
- END DO
- jm = nlay_s + nlay_i + 1
- DO ji = 1, npti
+ jm = nlay_s + nlay_i + 1
! ice bottom term
ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1)
ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
ztrid (ji,jm,3) = 0._wp
zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * &
& ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )
- END DO
- DO ji = 1, npti
! !---------------------!
IF( h_s_1d(ji) > 0._wp ) THEN ! snow-covered cells !
! !---------------------!
@@ -518,20 +492,19 @@ CONTAINS
zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &
& + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp
ENDIF
-
+
ENDIF
ENDIF
!
zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
zdiagbis(ji,jm_min(ji)) = ztrid (ji,jm_min(ji),2)
!
- END DO
- !
- !------------------------------
- ! 8) tridiagonal system solving
- !------------------------------
- ! Solve the tridiagonal system with Gauss elimination method.
- ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
+ !
+ !------------------------------
+ ! 8) tridiagonal system solving
+ !------------------------------
+ ! Solve the tridiagonal system with Gauss elimination method.
+ ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
!!$ jm_maxt = 0
!!$ jm_mint = nlay_i+5
!!$ DO ji = 1, npti
@@ -547,50 +520,39 @@ CONTAINS
!!$ zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1)
!!$ END DO
!!$ END DO
- ! clem: maybe one should find a way to reverse this loop for mpi performance
- DO ji = 1, npti
+ ! clem: maybe one should find a way to reverse this loop for mpi performance
jm_mint = jm_min(ji)
jm_maxt = jm_max(ji)
DO jm = jm_mint+1, jm_maxt
zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1)
zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1)
END DO
- END DO
- ! ice temperatures
- DO ji = 1, npti
+ ! ice temperatures
! Variable used after iterations
! Value must be frozen after convergence for MPP independance reason
IF ( .NOT. l_T_converged(ji) ) &
t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
- END DO
-
- DO jm = nlay_i + nlay_s, nlay_s + 2, -1
- DO ji = 1, npti
+
+ DO jm = nlay_i + nlay_s, nlay_s + 2, -1
jk = jm - nlay_s - 1
IF ( .NOT. l_T_converged(ji) ) &
t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
END DO
- END DO
- ! snow temperatures
- DO ji = 1, npti
+ ! snow temperatures
! Variables used after iterations
! Value must be frozen after convergence for MPP independance reason
IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
& t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
- END DO
- !!clem SNWLAY
- DO jm = nlay_s, 2, -1
- DO ji = 1, npti
+ !!clem SNWLAY
+ DO jm = nlay_s, 2, -1
jk = jm - 1
IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
& t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm)
END DO
- END DO
- ! surface temperature
- DO ji = 1, npti
+ ! surface temperature
IF( .NOT. l_T_converged(ji) ) THEN
ztsub(ji) = t_su_1d(ji)
IF( t_su_1d(ji) < rt0 ) THEN
@@ -598,20 +560,17 @@ CONTAINS
& ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
ENDIF
ENDIF
- END DO
- !
- !--------------------------------------------------------------
- ! 9) Has the scheme converged?, end of the iterative procedure
- !--------------------------------------------------------------
- ! check that nowhere it has started to melt
- ! zdti_max is a measure of error, it has to be under zdti_bnd
-
- DO ji = 1, npti
-
+ !
+ !--------------------------------------------------------------
+ ! 9) Has the scheme converged?, end of the iterative procedure
+ !--------------------------------------------------------------
+ ! check that nowhere it has started to melt
+ ! zdti_max is a measure of error, it has to be under zdti_bnd
+
zdti_max = 0._wp
IF ( .NOT. l_T_converged(ji) ) THEN
-
+
t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )
@@ -633,58 +592,53 @@ CONTAINS
tice_cvgerr_1d(ji) = zdti_max
tice_cvgstp_1d(ji) = REAL(iconv)
ENDIF
-
+
IF( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE.
-
+
ENDIF
-
- END DO
-
- !----------------------------------------!
- ! !
- ! Conduction flux is on !
- ! !
- !----------------------------------------!
- !
+
+ END DO
+ !
ELSEIF( k_cnd == np_cnd_ON ) THEN
+ !----------------------------------------!
+ ! !
+ ! Conduction flux is on !
+ ! !
+ !----------------------------------------!
!
! ==> we use a modified BL99 solver with conduction flux (qcn_ice) as forcing term
!
- !----------------------------
- ! 7) tridiagonal system terms
- !----------------------------
- ! layer denotes the number of the layer in the snow or in the ice
- ! jm denotes the reference number of the equation in the tridiagonal
- ! system, terms of tridiagonal system are indexed as following :
- ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
-
- ! ice interior terms (top equation has the same form as the others)
- ztrid (1:npti,:,:) = 0._wp
- zindterm(1:npti,:) = 0._wp
- zindtbis(1:npti,:) = 0._wp
- zdiagbis(1:npti,:) = 0._wp
-
- DO jm = nlay_s + 2, nlay_s + nlay_i
- DO ji = 1, npti
+ DO ji = 1, npti
+ !----------------------------
+ ! 7) tridiagonal system terms
+ !----------------------------
+ ! layer denotes the number of the layer in the snow or in the ice
+ ! jm denotes the reference number of the equation in the tridiagonal
+ ! system, terms of tridiagonal system are indexed as following :
+ ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
+
+ ! ice interior terms (top equation has the same form as the others)
+ ztrid (ji,:,:) = 0._wp
+ zindterm(ji,:) = 0._wp
+ zindtbis(ji,:) = 0._wp
+ zdiagbis(ji,:) = 0._wp
+
+ DO jm = nlay_s + 2, nlay_s + nlay_i
jk = jm - nlay_s - 1
ztrid (ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
ztrid (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
ztrid (ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk)
zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
END DO
- ENDDO
- jm = nlay_s + nlay_i + 1
- DO ji = 1, npti
! ice bottom term
+ jm = nlay_s + nlay_i + 1
ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1)
ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
ztrid (ji,jm,3) = 0._wp
zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * &
& ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )
- ENDDO
- DO ji = 1, npti
! !---------------------!
IF( h_s_1d(ji) > 0._wp ) THEN ! snow-covered cells !
! !---------------------!
@@ -738,13 +692,12 @@ CONTAINS
zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
zdiagbis(ji,jm_min(ji)) = ztrid (ji,jm_min(ji),2)
!
- END DO
- !
- !------------------------------
- ! 8) tridiagonal system solving
- !------------------------------
- ! Solve the tridiagonal system with Gauss elimination method.
- ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
+ !
+ !------------------------------
+ ! 8) tridiagonal system solving
+ !------------------------------
+ ! Solve the tridiagonal system with Gauss elimination method.
+ ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
!!$ jm_maxt = 0
!!$ jm_mint = nlay_i+5
!!$ DO ji = 1, npti
@@ -759,61 +712,49 @@ CONTAINS
!!$ zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
!!$ END DO
!!$ END DO
- ! clem: maybe one should find a way to reverse this loop for mpi performance
- DO ji = 1, npti
+ ! clem: maybe one should find a way to reverse this loop for mpi performance
jm_mint = jm_min(ji)
jm_maxt = jm_max(ji)
DO jm = jm_mint+1, jm_maxt
zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1)
zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1)
END DO
- END DO
- ! ice temperatures
- DO ji = 1, npti
+ ! ice temperatures
! Variable used after iterations
! Value must be frozen after convergence for MPP independance reason
IF ( .NOT. l_T_converged(ji) ) &
t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
- END DO
- DO jm = nlay_i + nlay_s, nlay_s + 2, -1
- DO ji = 1, npti
+ DO jm = nlay_i + nlay_s, nlay_s + 2, -1
IF ( .NOT. l_T_converged(ji) ) THEN
jk = jm - nlay_s - 1
t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
ENDIF
END DO
- END DO
- ! snow temperatures
- DO ji = 1, npti
+ ! snow temperatures
! Variables used after iterations
! Value must be frozen after convergence for MPP independance reason
IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
& t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
- END DO
- !!clem SNWLAY
- DO jm = nlay_s, 2, -1
- DO ji = 1, npti
+ !!clem SNWLAY
+ DO jm = nlay_s, 2, -1
jk = jm - 1
IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
& t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm)
END DO
- END DO
- !
- !--------------------------------------------------------------
- ! 9) Has the scheme converged?, end of the iterative procedure
- !--------------------------------------------------------------
- ! check that nowhere it has started to melt
- ! zdti_max is a measure of error, it has to be under zdti_bnd
-
- DO ji = 1, npti
-
+ !
+ !--------------------------------------------------------------
+ ! 9) Has the scheme converged?, end of the iterative procedure
+ !--------------------------------------------------------------
+ ! check that nowhere it has started to melt
+ ! zdti_max is a measure of error, it has to be under zdti_bnd
+
zdti_max = 0._wp
-
+
IF ( .NOT. l_T_converged(ji) ) THEN
-
+
IF( h_s_1d(ji) > 0._wp ) THEN
DO jk = 1, nlay_s
t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
--
GitLab