`_.
+| Start the update, save your changes and verify instantly the HTML export in your browser.
+
+.. warning::
+
+ | Your modifications are not taken into account?
+ | For symlink file, you will have to close it to update the HTML export. Otherwise look at the log of the Sphinx build, you probably made a typo!
+
+.. hint::
+
+ Are there broken links? Fix "Page not found" errors by running `make linkcheck`
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/doc/rst/requirements.txt
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/doc/rst/requirements.txt (revision 10901)
+++ (revision )
@@ -1,52 +1,0 @@
-alabaster==0.7.12
-Babel==2.6.0
-Brlapi==0.6.5
-certifi==2018.11.29
-chardet==3.0.4
-cryptography==1.7.1
-cupshelpers==1.0
-docutils==0.14
-httplib2==0.9.2
-idna==2.8
-imagesize==1.1.0
-Jinja2==2.10
-keyring==10.1
-keyrings.alt==1.3
-latexcodec==1.0.5
-louis==3.0.0
-Mako==1.0.6
-MarkupSafe==1.1.0
-oset==0.1.3
-packaging==18.0
-pexpect==4.2.1
-Pillow==4.0.0
-ptyprocess==0.5.1
-pyasn1==0.1.9
-pybtex==0.22.2
-pybtex-docutils==0.2.1
-pycrypto==2.6.1
-pycups==1.9.73
-pycurl==7.43.0
-Pygments==2.3.1
-pygobject==3.22.0
-pyparsing==2.3.1
-pysmbc==1.0.15.6
-python-apt==1.4.0b3
-python-debian==0.1.30
-python-debianbts==2.6.1
-pytz==2018.9
-pyxdg==0.25
-PyYAML==3.13
-reportbug==7.1.7
-reportlab==3.3.0
-requests==2.21.0
-roman==2.0.0
-SecretStorage==2.3.1
-six==1.12.0
-snowballstemmer==1.2.1
-Sphinx==1.8.3
-sphinx-rtd-theme==0.4.2
-sphinxcontrib-bibtex==0.4.2
-sphinxcontrib-websupport==1.1.0
-unattended-upgrades==0.1
-urllib3==1.24.1
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/doc/rst/source/_templates/layout.html
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/doc/rst/source/_templates/layout.html (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/doc/rst/source/_templates/layout.html (revision 11083)
@@ -12,5 +12,5 @@
Community ocean model for multifarious space and time scales
- {{ release }}
+ {{ version }}
{% include "searchbox.html" %}
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/doc/rst/source/conf.py
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/doc/rst/source/conf.py (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/doc/rst/source/conf.py (revision 11083)
@@ -92,4 +92,5 @@
# -- Customisation -----------------------------------------------------------
+# Timestamping
import datetime
year = datetime.date.today().year
@@ -108,2 +109,7 @@
# Include common directives for every rst file
rst_epilog = open('global.rst', 'r').read()
+
+# SVN revision
+import subprocess
+revision = subprocess.check_output("svnversion").decode("utf-8")
+rst_prolog = '.. |revision| replace:: %s' % revision
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icecor.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icecor.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icecor.F90 (revision 11083)
@@ -66,7 +66,6 @@
WRITE(numout,*) '~~~~~~~'
ENDIF
- !
! !-----------------------------------------------------
- ! ! ice thickness must exceed himin (for ice diff) !
+ ! ! ice thickness must exceed himin (for temp. diff.) !
! !-----------------------------------------------------
WHERE( a_i(:,:,:) >= epsi20 ) ; h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:)
@@ -97,5 +96,4 @@
END DO
ENDIF
-
! !-----------------------------------------------------
! ! Rebin categories with thickness out of bounds !
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icectl.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icectl.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icectl.F90 (revision 11083)
@@ -69,5 +69,5 @@
!!
REAL(wp) :: zv, zs, zt, zfs, zfv, zft
- REAL(wp) :: zvmin, zamin, zamax
+ REAL(wp) :: zvmin, zamin, zamax, zeimin, zesmin, zsmin
REAL(wp) :: zvtrp, zetrp
REAL(wp) :: zarea, zv_sill, zs_sill, zt_sill
@@ -141,7 +141,10 @@
zetrp = glob_sum( 'icectl', ( diag_trp_ei + diag_trp_es ) * e1e2t ) * zconv
- zvmin = glob_min( 'icectl', v_i )
- zamax = glob_max( 'icectl', SUM( a_i, dim=3 ) )
- zamin = glob_min( 'icectl', a_i )
+ zamax = glob_max( 'icectl', SUM( a_i, dim=3 ) )
+ zvmin = glob_min( 'icectl', v_i )
+ zamin = glob_min( 'icectl', a_i )
+ zsmin = glob_min( 'icectl', sv_i )
+ zeimin = glob_min( 'icectl', SUM( e_i, dim=3 ) )
+ zesmin = glob_min( 'icectl', SUM( e_s, dim=3 ) )
! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)
@@ -152,13 +155,18 @@
IF(lwp) THEN
- IF ( ABS( zv ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day] (',cd_routine,') = ',zv
- IF ( ABS( zs ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs
- IF ( ABS( zt ) > zt_sill ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',zt
- IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',zvmin
- IF ( zamax > MAX( rn_amax_n, rn_amax_s ) + epsi10 &
- & .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' .AND. cd_routine /= 'Hbig' ) &
- & WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax
- IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin
-!clem: the following check fails when using UMx advection scheme (see comments in icedyn_adv.F90)
+ ! check conservation issues
+ IF ( ABS( zv ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day] (',cd_routine,') = ',zv
+ IF ( ABS( zs ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs
+ IF ( ABS( zt ) > zt_sill ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',zt
+ ! check maximum ice concentration
+ IF ( zamax > MAX( rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) &
+ & WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax
+ ! check negative values
+ IF ( zvmin < 0. ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',zvmin
+ IF ( zamin < 0. ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin
+ IF ( zsmin < 0. ) WRITE(numout,*) 'violation s_i<0 (',cd_routine,') = ',zsmin
+ IF ( zeimin < 0. ) WRITE(numout,*) 'violation e_i<0 (',cd_routine,') = ',zeimin
+ IF ( zesmin < 0. ) WRITE(numout,*) 'violation e_s<0 (',cd_routine,') = ',zesmin
+!clem: the following check fails (I think...)
! IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN
! WRITE(numout,*) 'violation vtrp [Mt/day] (',cd_routine,') = ',zvtrp
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn.F90 (revision 11083)
@@ -74,8 +74,7 @@
INTEGER, INTENT(in) :: kt ! ice time step
!!
- INTEGER :: ji, jj, jl ! dummy loop indices
+ INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zcoefu, zcoefv
- REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max
- REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i
!!--------------------------------------------------------------------
!
@@ -89,48 +88,38 @@
ENDIF
!
- IF( ln_landfast_home ) THEN !-- Landfast ice parameterization
- tau_icebfr(:,:) = 0._wp
- DO jl = 1, jpl
- WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
- END DO
- ENDIF
- !
- ! !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig)
- DO jl = 1, jpl
- DO jj = 2, jpjm1
- DO ji = fs_2, fs_jpim1
- zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj ,jl), h_ip_b(ji ,jj+1,jl), &
- & h_ip_b(ji-1,jj ,jl), h_ip_b(ji ,jj-1,jl), &
- & h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), &
- & h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) )
- zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj ,jl), h_i_b (ji ,jj+1,jl), &
- & h_i_b (ji-1,jj ,jl), h_i_b (ji ,jj-1,jl), &
- & h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), &
- & h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) )
- zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj ,jl), h_s_b (ji ,jj+1,jl), &
- & h_s_b (ji-1,jj ,jl), h_s_b (ji ,jj-1,jl), &
- & h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), &
- & h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) )
- END DO
- END DO
- END DO
- CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. )
- !
- !
- SELECT CASE( nice_dyn ) !-- Set which dynamics is running
+ ! retrieve thickness from volume for landfast param. and UMx advection scheme
+ WHERE( a_i(:,:,:) >= epsi20 )
+ h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:)
+ h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:)
+ ELSEWHERE
+ h_i(:,:,:) = 0._wp
+ h_s(:,:,:) = 0._wp
+ END WHERE
+ !
+ WHERE( a_ip(:,:,:) >= epsi20 )
+ h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)
+ ELSEWHERE
+ h_ip(:,:,:) = 0._wp
+ END WHERE
+ !
+ !
+ SELECT CASE( nice_dyn ) !-- Set which dynamics is running
CASE ( np_dynALL ) !== all dynamical processes ==!
- CALL ice_dyn_rhg ( kt ) ! -- rheology
- CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness
- CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting
- CALL ice_cor ( kt , 1 ) ! -- Corrections
-
+ !
+ CALL ice_dyn_rhg ( kt ) ! -- rheology
+ CALL ice_dyn_adv ( kt ) ! -- advection of ice
+ CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting
+ CALL ice_cor ( kt , 1 ) ! -- Corrections
+ !
CASE ( np_dynRHGADV ) !== no ridge/raft & no corrections ==!
- CALL ice_dyn_rhg ( kt ) ! -- rheology
- CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness
- CALL Hpiling ! -- simple pile-up (replaces ridging/rafting)
-
+ !
+ CALL ice_dyn_rhg ( kt ) ! -- rheology
+ CALL ice_dyn_adv ( kt ) ! -- advection of ice
+ CALL Hpiling ! -- simple pile-up (replaces ridging/rafting)
+ CALL ice_var_zapsmall ! -- zap small areas
+ !
CASE ( np_dynADV1D ) !== pure advection ==! (1D)
- ALLOCATE( zdivu_i(jpi,jpj) )
+ !
! --- monotonicity test from Schar & Smolarkiewicz 1996 --- !
! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
@@ -145,20 +134,8 @@
END DO
! ---
- CALL ice_dyn_adv ( kt ) ! -- advection of ice
-
- ! diagnostics: divergence at T points
- DO jj = 2, jpjm1
- DO ji = 2, jpim1
- zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
- & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
- END DO
- END DO
- CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
- IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:) )
-
- DEALLOCATE( zdivu_i )
-
+ CALL ice_dyn_adv ( kt ) ! -- advection of ice
+ !
CASE ( np_dynADV2D ) !== pure advection ==! (2D w prescribed velocities)
- ALLOCATE( zdivu_i(jpi,jpj) )
+ !
u_ice(:,:) = rn_uice * umask(:,:,1)
v_ice(:,:) = rn_vice * vmask(:,:,1)
@@ -166,19 +143,30 @@
!CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
! ---
- CALL ice_dyn_adv ( kt ) ! -- advection of ice
-
- ! diagnostics: divergence at T points
- DO jj = 2, jpjm1
- DO ji = 2, jpim1
- zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
- & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
+ CALL ice_dyn_adv ( kt ) ! -- advection of ice
+
+ END SELECT
+ !
+ !
+ ! diagnostics: divergence at T points
+ IF( iom_use('icediv') ) THEN
+ !
+ SELECT CASE( nice_dyn )
+
+ CASE ( np_dynADV1D , np_dynADV2D )
+
+ ALLOCATE( zdivu_i(jpi,jpj) )
+ DO jj = 2, jpjm1
+ DO ji = 2, jpim1
+ zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
+ & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
+ END DO
END DO
- END DO
- CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
- IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:) )
-
- DEALLOCATE( zdivu_i )
-
- END SELECT
+ CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
+ CALL iom_put( "icediv" , zdivu_i(:,:) )
+ DEALLOCATE( zdivu_i )
+
+ END SELECT
+ !
+ ENDIF
!
! controls
@@ -188,92 +176,4 @@
- SUBROUTINE Hbig( phi_max, phs_max, phip_max )
- !!-------------------------------------------------------------------
- !! *** ROUTINE Hbig ***
- !!
- !! ** Purpose : Thickness correction in case advection scheme creates
- !! abnormally tick ice or snow
- !!
- !! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points
- !! (before advection) and reduce it by adapting ice concentration
- !! 2- check whether snow thickness is larger than the surrounding 9-points
- !! (before advection) and reduce it by sending the excess in the ocean
- !! 3- check whether snow load deplets the snow-ice interface below sea level$
- !! and reduce it by sending the excess in the ocean
- !! 4- correct pond fraction to avoid a_ip > a_i
- !!
- !! ** input : Max thickness of the surrounding 9-points
- !!-------------------------------------------------------------------
- REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts
- !
- INTEGER :: ji, jj, jl ! dummy loop indices
- REAL(wp) :: zhip, zhi, zhs, zvs_excess, zfra
- !!-------------------------------------------------------------------
- ! controls
- IF( ln_icediachk ) CALL ice_cons_hsm(0, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
- !
- CALL ice_var_zapsmall !-- zap small areas
- !
- DO jl = 1, jpl
- DO jj = 1, jpj
- DO ji = 1, jpi
- IF ( v_i(ji,jj,jl) > 0._wp ) THEN
- !
- ! ! -- check h_ip -- !
- ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
- IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN
- zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) )
- IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN
- a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl)
- ENDIF
- ENDIF
- !
- ! ! -- check h_i -- !
- ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
- zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
- IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
- a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m)
- ENDIF
- !
- ! ! -- check h_s -- !
- ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
- zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
- IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
- zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
- !
- wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice
- hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
- !
- e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
- v_s(ji,jj,jl) = a_i(ji,jj,jl) * phs_max(ji,jj,jl)
- ENDIF
- !
- ! ! -- check snow load -- !
- ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
- ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
- ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
- zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
- IF( zvs_excess > 0._wp ) THEN
- zfra = ( v_s(ji,jj,jl) - zvs_excess ) / MAX( v_s(ji,jj,jl), epsi20 )
- wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice
- hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
- !
- e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
- v_s(ji,jj,jl) = v_s(ji,jj,jl) - zvs_excess
- ENDIF
-
- ENDIF
- END DO
- END DO
- END DO
- ! !-- correct pond fraction to avoid a_ip > a_i
- WHERE( a_ip(:,:,:) > a_i(:,:,:) ) a_ip(:,:,:) = a_i(:,:,:)
- !
- ! controls
- IF( ln_icediachk ) CALL ice_cons_hsm(1, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
- !
- END SUBROUTINE Hbig
-
-
SUBROUTINE Hpiling
!!-------------------------------------------------------------------
@@ -290,6 +190,4 @@
! controls
IF( ln_icediachk ) CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
- !
- CALL ice_var_zapsmall !-- zap small areas
!
at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
@@ -337,16 +235,16 @@
WRITE(numout,*) '~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namdyn:'
- WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) ln_dynALL = ', ln_dynALL
- WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) ln_dynRHGADV = ', ln_dynRHGADV
- WRITE(numout,*) ' Advection 1D only (Schar & Smolarkiewicz 1996) ln_dynADV1D = ', ln_dynADV1D
- WRITE(numout,*) ' Advection 2D only (rn_uvice + adv) ln_dynADV2D = ', ln_dynADV2D
- WRITE(numout,*) ' with prescribed velocity given by (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')'
- WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat
- WRITE(numout,*) ' Landfast: param from Lemieux 2016 ln_landfast_L16 = ', ln_landfast_L16
- WRITE(numout,*) ' Landfast: param from home made ln_landfast_home= ', ln_landfast_home
- WRITE(numout,*) ' fraction of ocean depth that ice must reach rn_depfra = ', rn_depfra
- WRITE(numout,*) ' maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr
- WRITE(numout,*) ' relax time scale (s-1) to reach static friction rn_lfrelax = ', rn_lfrelax
- WRITE(numout,*) ' isotropic tensile strength rn_tensile = ', rn_tensile
+ WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) ln_dynALL = ', ln_dynALL
+ WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) ln_dynRHGADV = ', ln_dynRHGADV
+ WRITE(numout,*) ' Advection 1D only (Schar & Smolarkiewicz 1996) ln_dynADV1D = ', ln_dynADV1D
+ WRITE(numout,*) ' Advection 2D only (rn_uvice + adv) ln_dynADV2D = ', ln_dynADV2D
+ WRITE(numout,*) ' with prescribed velocity given by (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',',rn_vice,')'
+ WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat
+ WRITE(numout,*) ' Landfast: param from Lemieux 2016 ln_landfast_L16 = ', ln_landfast_L16
+ WRITE(numout,*) ' Landfast: param from home made ln_landfast_home= ', ln_landfast_home
+ WRITE(numout,*) ' fraction of ocean depth that ice must reach rn_depfra = ', rn_depfra
+ WRITE(numout,*) ' maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr
+ WRITE(numout,*) ' relax time scale (s-1) to reach static friction rn_lfrelax = ', rn_lfrelax
+ WRITE(numout,*) ' isotropic tensile strength rn_tensile = ', rn_tensile
WRITE(numout,*)
ENDIF
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_adv.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_adv.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_adv.F90 (revision 11083)
@@ -64,10 +64,9 @@
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! number of iteration
- !
- INTEGER :: jl ! dummy loop indice
- REAL(wp), DIMENSION(jpi,jpj) :: zmask ! fraction of time step with v_i < 0
!!---------------------------------------------------------------------
!
- IF( ln_timing ) CALL timing_start('icedyn_adv')
+ ! controls
+ IF( ln_timing ) CALL timing_start('icedyn_adv') ! timing
+ IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
!
IF( kt == nit000 .AND. lwp ) THEN
@@ -76,44 +75,20 @@
WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
-
- ! conservation test
- IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
-
- !----------
- ! Advection
- !----------
+ !
+ !---------------!
+ !== Advection ==!
+ !---------------!
SELECT CASE( nice_adv )
! !-----------------------!
CASE( np_advUMx ) ! ULTIMATE-MACHO scheme !
! !-----------------------!
- CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
- ! !-----------------------!
+ CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, &
+ & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+ ! !-----------------------!
CASE( np_advPRA ) ! PRATHER scheme !
! !-----------------------!
- CALL ice_dyn_adv_pra( kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+ CALL ice_dyn_adv_pra( kt, u_ice, v_ice, &
+ & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
END SELECT
-
- !----------------------------
- ! Debug the advection schemes
- !----------------------------
- ! clem: At least one advection scheme above is not strictly positive => UMx
- ! In Prather, I am not sure if the fields are bounded by 0 or not (it seems yes)
- ! In UMx , advected fields are not perfectly bounded and negative values can appear.
- ! These values are usually very small but in some occasions they can also be non-negligible
- ! Therefore one needs to bound the advected fields by 0 (though this is not a clean fix)
- !
- ! record the negative values resulting from UMx
- zmask(:,:) = 0._wp ! keep the init to 0 here
- DO jl = 1, jpl
- WHERE( v_i(:,:,jl) < 0._wp ) zmask(:,:) = 1._wp
- END DO
- IF( iom_use('iceneg_pres') ) CALL iom_put("iceneg_pres", zmask ) ! fraction of time step with v_i < 0
- IF( iom_use('iceneg_volu') ) CALL iom_put("iceneg_volu", SUM(MIN( v_i, 0. ), dim=3 ) ) ! negative ice volume (only)
- IF( iom_use('iceneg_hfx' ) ) CALL iom_put("iceneg_hfx" , ( SUM(SUM( MIN( e_i(:,:,1:nlay_i,:), 0. ) & ! negative ice heat content (only)
- & , dim=4 ), dim=3 ) )* r1_rdtice ) ! -- eq. heat flux [W/m2]
- !
- ! ==> conservation is ensured by calling this routine below,
- ! however the global ice volume is then changed by advection (but errors are small)
- CALL ice_var_zapneg( ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
!------------
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_adv_umx.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_adv_umx.F90 (revision 11083)
@@ -11,8 +11,8 @@
!! 'key_si3' SI3 sea-ice model
!!----------------------------------------------------------------------
- !! ice_dyn_adv_umx : update the tracer trend with the 3D advection trends using a TVD scheme
+ !! ice_dyn_adv_umx : update the tracer fields
!! ultimate_x(_y) : compute a tracer value at velocity points using ULTIMATE scheme at various orders
- !! macho : ???
- !! nonosc_ice : compute monotonic tracer fluxes by a non-oscillatory algorithm
+ !! macho : compute the fluxes
+ !! nonosc_ice : limit the fluxes using a non-oscillatory algorithm
!!----------------------------------------------------------------------
USE phycst ! physical constant
@@ -32,27 +32,22 @@
PUBLIC ice_dyn_adv_umx ! called by icedyn_adv.F90
-
- REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6
- REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120
-
- ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid)
- INTEGER :: kn_limiter = 1
-
- ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme
- ! clem: if set to true, the 2D test case "diagonal advection" does not work (I do not understand why)
- ! but in realistic cases, it avoids having very negative ice temperature (-50) at low ice concentration
- LOGICAL :: ll_neg = .TRUE.
-
- ! alternate directions for upstream
- LOGICAL :: ll_upsxy = .TRUE.
-
- ! alternate directions for high order
- LOGICAL :: ll_hoxy = .TRUE.
-
- ! prelimiter: use it to avoid overshoot in H
- ! clem: if set to true, the 2D test case "diagnoal advection" does not work (I do not understand why)
- LOGICAL :: ll_prelimiter_zalesak = .FALSE. ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D
-
-
+ !
+ INTEGER, PARAMETER :: np_advS = 1 ! advection for S and T: dVS/dt = -div( uVS ) => np_advS = 1
+ ! or dVS/dt = -div( uA * uHS / u ) => np_advS = 2
+ ! or dVS/dt = -div( uV * uS / u ) => np_advS = 3
+ INTEGER, PARAMETER :: np_limiter = 1 ! limiter: 1 = nonosc
+ ! 2 = superbee
+ ! 3 = h3
+ LOGICAL :: ll_upsxy = .TRUE. ! alternate directions for upstream
+ LOGICAL :: ll_hoxy = .TRUE. ! alternate directions for high order
+ LOGICAL :: ll_neg = .TRUE. ! if T interpolated at u/v points is negative or v_i < 1.e-6
+ ! then interpolate T at u/v points using the upstream scheme
+ LOGICAL :: ll_prelim = .FALSE. ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D
+ !
+ REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6
+ REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120
+ !
+ INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: imsk_small, jmsk_small
+ !
!! * Substitutions
# include "vectopt_loop_substitute.h90"
@@ -64,5 +59,5 @@
CONTAINS
- SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, &
+ SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, &
& pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
!!----------------------------------------------------------------------
@@ -79,4 +74,7 @@
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity
+ REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_i ! ice thickness
+ REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_s ! snw thickness
+ REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_ip ! ice pond thickness
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
@@ -93,18 +91,43 @@
INTEGER :: icycle ! number of sub-timestep for the advection
REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers
- REAL(wp) :: zdt
- REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! send zcflnow and receive zcflprv
- REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box
+ REAL(wp) :: zdt, zvi_cen
+ REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication
+ REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box
REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2
- REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho
- REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai, z1_aip
- REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhvar
+ REAL(wp), DIMENSION(jpi,jpj,jpl) :: zu_cat, zv_cat
+ REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho, zua_ups, zva_ups
+ REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai , z1_aip, zhvar
+ REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max
+ !
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs
!!----------------------------------------------------------------------
!
IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
!
- ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !
- ! When needed, the advection split is applied at the next time-step in order to avoid blocking global comm.
- ! ...this should not affect too much the stability... Was ok on the tests we did...
+ ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- !
+ DO jl = 1, jpl
+ DO jj = 2, jpjm1
+ DO ji = fs_2, fs_jpim1
+ zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), &
+ & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), &
+ & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
+ & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
+ zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), &
+ & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), &
+ & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
+ & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
+ zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), &
+ & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), &
+ & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
+ & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
+ END DO
+ END DO
+ END DO
+ CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. )
+ !
+ !
+ ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
+ ! Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
+ ! this should not affect too much the stability
zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
@@ -116,5 +139,4 @@
ELSE ; icycle = 1
ENDIF
-
zdt = rdt_ice / REAL(icycle)
@@ -122,5 +144,11 @@
zudy(:,:) = pu_ice(:,:) * e2u(:,:)
zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
-
+ !
+ ! setup transport for each ice cat
+ DO jl = 1, jpl
+ zu_cat(:,:,jl) = zudy(:,:)
+ zv_cat(:,:,jl) = zvdx(:,:)
+ END DO
+ !
! --- define velocity for advection: u*grad(H) --- !
DO jj = 2, jpjm1
@@ -154,73 +182,180 @@
END WHERE
!
- ! set u*a=u for advection of A only
- DO jl = 1, jpl
- zua_ho(:,:,jl) = zudy(:,:)
- zva_ho(:,:,jl) = zvdx(:,:)
- END DO
-
+ ! setup a mask where advection will be upstream
+ IF( ll_neg ) THEN
+ IF( .NOT. ALLOCATED(imsk_small) ) ALLOCATE( imsk_small(jpi,jpj,jpl) )
+ IF( .NOT. ALLOCATED(jmsk_small) ) ALLOCATE( jmsk_small(jpi,jpj,jpl) )
+ DO jl = 1, jpl
+ DO jj = 1, jpjm1
+ DO ji = 1, jpim1
+ zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) )
+ IF( zvi_cen < epsi06) THEN ; imsk_small(ji,jj,jl) = 0
+ ELSE ; imsk_small(ji,jj,jl) = 1 ; ENDIF
+ zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) )
+ IF( zvi_cen < epsi06) THEN ; jmsk_small(ji,jj,jl) = 0
+ ELSE ; jmsk_small(ji,jj,jl) = 1 ; ENDIF
+ END DO
+ END DO
+ END DO
+ ENDIF
+ !
+ ! ----------------------- !
+ ! ==> start advection <== !
+ ! ----------------------- !
+ !
+ !== Ice area ==!
zamsk = 1._wp
- CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) !-- Ice area
- zamsk = 0._wp
- !
- zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
- CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i ) !-- Ice volume
- !
- zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
- CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s ) !-- Snw volume
- !
- zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
- CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i ) !-- Salt content
- !
- DO jk = 1, nlay_i
- zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
- CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) ) !-- Ice heat content
- END DO
- !
- DO jk = 1, nlay_s
- zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
- CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) ) !-- Snw heat content
- END DO
- !
- IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN !-- Ice Age
- ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with
- ! fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that
- ! spread into the domain. Therefore we cheat and consider that ice age should be advected as ice concentration
- !!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
- !!CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i )
- ! set u*a=u for advection of ice age
- DO jl = 1, jpl
- zua_ho(:,:,jl) = zudy(:,:)
- zva_ho(:,:,jl) = zvdx(:,:)
- END DO
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zu_cat , zv_cat , zcu_box, zcv_box, &
+ & pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho )
+ !
+ ! ! --------------------------------- !
+ IF( np_advS == 1 ) THEN ! -- advection form: -div( uVS ) -- !
+ ! ! --------------------------------- !
+ zamsk = 0._wp
+ !== Ice volume ==!
+ zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+ & zhvar, pv_i, zua_ups, zva_ups )
+ !== Snw volume ==!
+ zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+ & zhvar, pv_s, zua_ups, zva_ups )
+ !
zamsk = 1._wp
- CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, poa_i, poa_i )
+ !== Salt content ==!
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+ & psv_i, psv_i )
+ !== Ice heat content ==!
+ DO jk = 1, nlay_i
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+ & pe_i(:,:,jk,:), pe_i(:,:,jk,:) )
+ END DO
+ !== Snw heat content ==!
+ DO jk = 1, nlay_s
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+ & pe_s(:,:,jk,:), pe_s(:,:,jk,:) )
+ END DO
+ !
+ ! ! ------------------------------------------ !
+ ELSEIF( np_advS == 2 ) THEN ! -- advection form: -div( uA * uHS / u ) -- !
+ ! ! ------------------------------------------ !
zamsk = 0._wp
+ !== Ice volume ==!
+ zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+ & zhvar, pv_i, zua_ups, zva_ups )
+ !== Snw volume ==!
+ zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+ & zhvar, pv_s, zua_ups, zva_ups )
+ !== Salt content ==!
+ zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
+ & zhvar, psv_i, zua_ups, zva_ups )
+ !== Ice heat content ==!
+ DO jk = 1, nlay_i
+ zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
+ & zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups )
+ END DO
+ !== Snw heat content ==!
+ DO jk = 1, nlay_s
+ zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
+ & zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups )
+ END DO
+ !
+ ! ! ----------------------------------------- !
+ ELSEIF( np_advS == 3 ) THEN ! -- advection form: -div( uV * uS / u ) -- !
+ ! ! ----------------------------------------- !
+ zamsk = 0._wp
+ !
+ ALLOCATE( zuv_ho (jpi,jpj,jpl), zvv_ho (jpi,jpj,jpl), &
+ & zuv_ups(jpi,jpj,jpl), zvv_ups(jpi,jpj,jpl), z1_vi(jpi,jpj,jpl), z1_vs(jpi,jpj,jpl) )
+ !
+ ! inverse of Vi
+ WHERE( pv_i(:,:,:) >= epsi20 ) ; z1_vi(:,:,:) = 1._wp / pv_i(:,:,:)
+ ELSEWHERE ; z1_vi(:,:,:) = 0.
+ END WHERE
+ ! inverse of Vs
+ WHERE( pv_s(:,:,:) >= epsi20 ) ; z1_vs(:,:,:) = 1._wp / pv_s(:,:,:)
+ ELSEWHERE ; z1_vs(:,:,:) = 0.
+ END WHERE
+ !
+ ! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays)
+ !
+ !== Ice volume ==!
+ zuv_ups = zua_ups
+ zvv_ups = zva_ups
+ zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+ & zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
+ !== Salt content ==!
+ zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, &
+ & zhvar, psv_i, zuv_ups, zvv_ups )
+ !== Ice heat content ==!
+ DO jk = 1, nlay_i
+ zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
+ & zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups )
+ END DO
+ !== Snow volume ==!
+ zuv_ups = zua_ups
+ zvv_ups = zva_ups
+ zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
+ & zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
+ !== Snw heat content ==!
+ DO jk = 1, nlay_s
+ zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:)
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
+ & zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups )
+ END DO
+ !
+ DEALLOCATE( zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs )
+ !
ENDIF
!
- IF ( ln_pnd_H12 ) THEN !-- melt ponds
- ! set u*a=u for advection of Ap only
- DO jl = 1, jpl
- zua_ho(:,:,jl) = zudy(:,:)
- zva_ho(:,:,jl) = zvdx(:,:)
- END DO
- !
+ !== Ice age ==!
+ IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN
zamsk = 1._wp
- CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! fraction
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+ & poa_i, poa_i )
+ ENDIF
+ !
+ !== melt ponds ==!
+ IF ( ln_pnd_H12 ) THEN
+ ! fraction
+ zamsk = 1._wp
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, &
+ & pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho )
+ ! volume
zamsk = 0._wp
- !
zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:)
- CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip ) ! volume
+ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
+ & zhvar, pv_ip, zua_ups, zva_ups )
ENDIF
!
+ !== Open water area ==!
zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
DO jj = 2, jpjm1
DO ji = fs_2, fs_jpim1
- pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !-- Open water area
+ pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &
& - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
END DO
END DO
- CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T', 1. )
- !
+ CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1. )
+ !
+ !
+ ! --- Ensure non-negative fields and in-bound thicknesses --- !
+ ! Remove negative values (conservation is ensured)
+ ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
+ CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+ !
+ ! Make sure ice thickness is not too big
+ ! (because ice thickness can be too large where ice concentration is very small)
+ CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+
END DO
!
@@ -228,5 +363,6 @@
- SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho )
+ SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, &
+ & pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho )
!!----------------------------------------------------------------------
!! *** ROUTINE adv_umx ***
@@ -235,40 +371,51 @@
!! tracers and add it to the general trend of tracer equations
!!
- !! ** Method : - calculate upstream fluxes and upstream solution for tracer H
+ !! ** Method : - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc
!! - calculate tracer H at u and v points (Ultimate)
- !! - calculate the high order fluxes using alterning directions (Macho?)
+ !! - calculate the high order fluxes using alterning directions (Macho)
!! - apply a limiter on the fluxes (nonosc_ice)
- !! - convert this tracer flux to a tracer content flux (uH -> uV)
- !! - calculate the high order solution for tracer content V
+ !! - convert this tracer flux to a "volume" flux (uH -> uV)
+ !! - apply a limiter a second time on the volumes fluxes (nonosc_ice)
+ !! - calculate the high order solution for V
!!
- !! ** Action : solve 2 equations => a) da/dt = -div(ua)
- !! b) dV/dt = -div(uV) using dH/dt = -u.grad(H)
- !! in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc_ice). This step is necessary to get a good H.
- !! - then we convert this flux to a "volume" flux this way => uH*ua/u
- !! where ua is the flux from eq. a)
- !! - at last we estimate dV/dt = -div(uH*ua/u)
+ !! ** Action : solve 3 equations => a) dA/dt = -div(uA)
+ !! b) dV/dt = -div(uV) using dH/dt = -u.grad(H)
+ !! c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S)
!!
- !! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc.
- !! - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now
- !! is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of
- !! the scheme automatically by applying a mask of the ice cover inside Ultimate (not done).
+ !! in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H.
+ !! - then we convert this flux to a "volume" flux this way => uH * uA / u
+ !! where uA is the flux from eq. a)
+ !! this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur)
+ !! - at last we estimate dV/dt = -div(uH * uA / u)
+ !!
+ !! in eq. c), one can solve the equation for S (ln_advS=T), then dVS/dt = -div(uV * uS / u)
+ !! or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u)
+ !!
+ !! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc.
+ !! - At the ice edge, Ultimate scheme can lead to:
+ !! 1) negative interpolated tracers at u-v points
+ !! 2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward
+ !! Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of
+ !! the scheme automatically by applying a mask of the ice cover inside Ultimate (not done).
+ !! Solution for 2): we set it to 0 in this case
!! - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good.
!! Large values of H can appear for very small ice concentration, and when it does it messes the things up since we
- !! work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D.
+ !! work on H (and not V). It is partly related to the multi-category approach
!! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice
- !! concentration is small).
- !! To-do: expand the prelimiter from zalesak to make it work in 2D
- !!----------------------------------------------------------------------
- REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
- INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
- INTEGER , INTENT(in ) :: jt ! number of sub-iteration
- INTEGER , INTENT(in ) :: kt ! number of iteration
- REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
- REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2
- REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u
- REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity
- REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field
- REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field
- REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes
+ !! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter
+ !! since sv_i and e_i are still good.
+ !!----------------------------------------------------------------------
+ REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
+ INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
+ INTEGER , INTENT(in ) :: jt ! number of sub-iteration
+ INTEGER , INTENT(in ) :: kt ! number of iteration
+ REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
+ REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2
+ REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u
+ REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity
+ REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field
+ REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field
+ REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL :: pua_ups, pva_ups ! upstream u*a fluxes
+ REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes
!
INTEGER :: ji, jj, jl ! dummy loop indices
@@ -303,7 +450,7 @@
DO jj = 1, jpjm1
DO ji = 1, fs_jpim1
- IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
- zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
- zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
+ IF( ABS( pu(ji,jj) ) > epsi10 ) THEN
+ zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj)
+ zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj)
ELSE
zfu_ho (ji,jj,jl) = 0._wp
@@ -311,7 +458,7 @@
ENDIF
!
- IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
- zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
- zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
+ IF( ABS( pv(ji,jj) ) > epsi10 ) THEN
+ zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj)
+ zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj)
ELSE
zfv_ho (ji,jj,jl) = 0._wp
@@ -321,15 +468,36 @@
END DO
END DO
+
+ ! the new "volume" fluxes must also be "flux corrected"
+ ! thus we calculate the upstream solution and apply a limiter again
+ DO jl = 1, jpl
+ DO jj = 2, jpjm1
+ DO ji = fs_2, fs_jpim1
+ ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) )
+ !
+ zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
+ END DO
+ END DO
+ END DO
+ CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. )
+ !
+ IF ( np_limiter == 1 ) THEN
+ CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
+ ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
+ CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho )
+ CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho )
+ ENDIF
+ !
ENDIF
- ! --ho
- ! in case of advection of A: output u*a
- ! -------------------------------------
+ ! --ho --ups
+ ! in case of advection of A: output u*a and u*a
+ ! -----------------------------------------------
IF( PRESENT( pua_ho ) ) THEN
DO jl = 1, jpl
DO jj = 1, jpjm1
DO ji = 1, fs_jpim1
- pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl)
- pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl)
- END DO
+ pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
+ pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
+ END DO
END DO
END DO
@@ -499,7 +667,7 @@
END DO
!
- IF ( kn_limiter == 1 ) THEN
+ IF ( np_limiter == 1 ) THEN
CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
- ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN
+ ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
@@ -517,5 +685,5 @@
END DO
END DO
- IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+ IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
DO jl = 1, jpl !-- first guess of tracer from u-flux
@@ -538,5 +706,5 @@
END DO
END DO
- IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+ IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
ELSE !== even ice time step: adv_y then adv_x ==!
@@ -549,5 +717,5 @@
END DO
END DO
- IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+ IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
!
DO jl = 1, jpl !-- first guess of tracer from v-flux
@@ -570,8 +738,8 @@
END DO
END DO
- IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+ IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
ENDIF
- IF( kn_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+ IF( np_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
ENDIF
@@ -609,7 +777,7 @@
!
! !-- ultimate interpolation of pt at u-point --!
- CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho )
+ CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho )
! !-- limiter in x --!
- IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+ IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
! !-- advective form update in zpt --!
DO jl = 1, jpl
@@ -619,5 +787,5 @@
& + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) &
& * pamsk &
- & ) * pdt ) * tmask(ji,jj,1)
+ & ) * pdt ) * tmask(ji,jj,1)
END DO
END DO
@@ -627,10 +795,10 @@
! !-- ultimate interpolation of pt at v-point --!
IF( ll_hoxy ) THEN
- CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
+ CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
ELSE
- CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho )
+ CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho )
ENDIF
! !-- limiter in y --!
- IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+ IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
!
!
@@ -638,7 +806,7 @@
!
! !-- ultimate interpolation of pt at v-point --!
- CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho )
+ CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho )
! !-- limiter in y --!
- IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
+ IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
! !-- advective form update in zpt --!
DO jl = 1, jpl
@@ -656,19 +824,19 @@
! !-- ultimate interpolation of pt at u-point --!
IF( ll_hoxy ) THEN
- CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
+ CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
ELSE
- CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho )
+ CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho )
ENDIF
! !-- limiter in x --!
- IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
+ IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
!
ENDIF
- IF( kn_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+ IF( np_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
!
END SUBROUTINE macho
- SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho )
+ SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho )
!!---------------------------------------------------------------------
!! *** ROUTINE ultimate_x ***
@@ -680,4 +848,5 @@
!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
!!----------------------------------------------------------------------
+ REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
@@ -806,5 +975,5 @@
DO jj = 1, jpjm1
DO ji = 1, fs_jpim1
- IF( pt_u(ji,jj,jl) < 0._wp ) THEN
+ IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
@@ -826,5 +995,5 @@
- SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho )
+ SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho )
!!---------------------------------------------------------------------
!! *** ROUTINE ultimate_y ***
@@ -836,4 +1005,5 @@
!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
!!----------------------------------------------------------------------
+ REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
@@ -959,5 +1129,5 @@
DO jj = 1, jpjm1
DO ji = 1, fs_jpim1
- IF( pt_v(ji,jj,jl) < 0._wp ) THEN
+ IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) &
& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
@@ -1023,5 +1193,5 @@
! | | | | *
! t_ups : i-1 i i+1 i+2
- IF( ll_prelimiter_zalesak ) THEN
+ IF( ll_prelim ) THEN
DO jl = 1, jpl
@@ -1102,10 +1272,11 @@
!
! ! up & down beta terms
- IF( zpos > 0._wp ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
- ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig
+ ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!)
+ IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
+ ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig
ENDIF
!
- IF( zneg > 0._wp ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
- ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig
+ IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
+ ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig
ENDIF
!
@@ -1149,18 +1320,4 @@
END DO
- ! clem test
-!! DO jj = 2, jpjm1
-!! DO ji = 2, fs_jpim1 ! vector opt.
-!! zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
-!! & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
-!! & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
-!! & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
-!! & ) * tmask(ji,jj,1)
-!! IF( zzt < -epsi20 ) THEN
-!! WRITE(numout,*) 'T<0 nonosc_ice',zzt
-!! ENDIF
-!! END DO
-!! END DO
-
END DO
!
@@ -1203,5 +1360,5 @@
Rjp = zslpx(ji+1,jj,jl)
- IF( kn_limiter == 3 ) THEN
+ IF( np_limiter == 3 ) THEN
IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm
@@ -1219,5 +1376,5 @@
pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
- ELSEIF( kn_limiter == 2 ) THEN
+ ELSEIF( np_limiter == 2 ) THEN
IF( Rj /= 0. ) THEN
IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj
@@ -1298,5 +1455,5 @@
Rjp = zslpy(ji,jj+1,jl)
- IF( kn_limiter == 3 ) THEN
+ IF( np_limiter == 3 ) THEN
IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm
@@ -1314,5 +1471,5 @@
pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
- ELSEIF( kn_limiter == 2 ) THEN
+ ELSEIF( np_limiter == 2 ) THEN
IF( Rj /= 0. ) THEN
@@ -1358,4 +1515,94 @@
END SUBROUTINE limiter_y
+
+ SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+ !!-------------------------------------------------------------------
+ !! *** ROUTINE Hbig ***
+ !!
+ !! ** Purpose : Thickness correction in case advection scheme creates
+ !! abnormally tick ice or snow
+ !!
+ !! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points
+ !! (before advection) and reduce it by adapting ice concentration
+ !! 2- check whether snow thickness is larger than the surrounding 9-points
+ !! (before advection) and reduce it by sending the excess in the ocean
+ !! 3- check whether snow load deplets the snow-ice interface below sea level$
+ !! and reduce it by sending the excess in the ocean
+ !! 4- correct pond fraction to avoid a_ip > a_i
+ !!
+ !! ** input : Max thickness of the surrounding 9-points
+ !!-------------------------------------------------------------------
+ REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
+ REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts
+ REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip
+ REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s
+ REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i
+ !
+ INTEGER :: ji, jj, jk, jl ! dummy loop indices
+ REAL(wp) :: z1_dt, zhip, zhi, zhs, zvs_excess, zfra
+ REAL(wp), DIMENSION(jpi,jpj) :: zswitch
+ !!-------------------------------------------------------------------
+ !
+ z1_dt = 1._wp / pdt
+ !
+ DO jl = 1, jpl
+
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+ !
+ ! ! -- check h_ip -- !
+ ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
+ IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
+ zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
+ IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
+ pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
+ ENDIF
+ ENDIF
+ !
+ ! ! -- check h_i -- !
+ ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
+ zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
+ IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+ pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m)
+ ENDIF
+ !
+ ! ! -- check h_s -- !
+ ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
+ zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
+ IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
+ zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
+ !
+ wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt
+ hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
+ !
+ pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
+ pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
+ ENDIF
+ !
+ ! ! -- check snow load -- !
+ ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
+ ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
+ ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
+ zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
+ IF( zvs_excess > 0._wp ) THEN
+ zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
+ wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
+ hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
+ !
+ pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
+ pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess
+ ENDIF
+
+ ENDIF
+ END DO
+ END DO
+ END DO
+ ! !-- correct pond fraction to avoid a_ip > a_i
+ WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:)
+ !
+ !
+ END SUBROUTINE Hbig
+
#else
!!----------------------------------------------------------------------
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rdgrft.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rdgrft.F90 (revision 11083)
@@ -142,8 +142,4 @@
INTEGER, PARAMETER :: jp_itermax = 20
!!-------------------------------------------------------------------
- ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem)
- ! likely due to truncation error ( i.e. 1. - 1. /= 0 )
- ! I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)
-
! controls
IF( ln_timing ) CALL timing_start('icedyn_rdgrft') ! timing
@@ -156,14 +152,14 @@
ENDIF
- CALL ice_var_zapsmall ! Zero out categories with very small areas
-
!--------------------------------
! 0) Identify grid cells with ice
!--------------------------------
+ at_i(:,:) = SUM( a_i, dim=3 )
+ !
npti = 0 ; nptidx(:) = 0
ipti = 0 ; iptidx(:) = 0
DO jj = 1, jpj
DO ji = 1, jpi
- IF ( at_i(ji,jj) > 0._wp ) THEN
+ IF ( at_i(ji,jj) > epsi10 ) THEN
npti = npti + 1
nptidx( npti ) = (jj - 1) * jpi + ji
@@ -178,10 +174,10 @@
! just needed here
- CALL tab_2d_1d( npti, nptidx(1:npti), zdivu(1:npti), divu_i(:,:) )
- CALL tab_2d_1d( npti, nptidx(1:npti), zdelt(1:npti), delta_i(:,:) )
+ CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i )
+ CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i )
! needed here and in the iteration loop
- 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_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i(:,:) )
+ 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_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i )
DO ji = 1, npti
@@ -310,6 +306,6 @@
! ! Ice thickness needed for rafting
- WHERE( pa_i(1:npti,:) > 0._wp ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
- ELSEWHERE ; zhi(1:npti,:) = 0._wp
+ WHERE( pa_i(1:npti,:) > epsi20 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
+ ELSEWHERE ; zhi(1:npti,:) = 0._wp
END WHERE
@@ -329,6 +325,6 @@
zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
!
- WHERE( zasum(1:npti) > 0._wp ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti)
- ELSEWHERE ; z1_asum(1:npti) = 0._wp
+ WHERE( zasum(1:npti) > epsi20 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti)
+ ELSEWHERE ; z1_asum(1:npti) = 0._wp
END WHERE
!
@@ -455,6 +451,6 @@
! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.
! NOTE: 0 < aksum <= 1
- WHERE( zaksum(1:npti) > 0._wp ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
- ELSEWHERE ; closing_gross(1:npti) = 0._wp
+ WHERE( zaksum(1:npti) > epsi20 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
+ ELSEWHERE ; closing_gross(1:npti) = 0._wp
END WHERE
@@ -466,5 +462,5 @@
DO ji = 1, npti
zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
- IF( zfac > pa_i(ji,jl) ) THEN
+ IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice
ENDIF
@@ -510,4 +506,5 @@
REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2
REAL(wp), DIMENSION(jpij) :: z1_ai ! 1 / a
+ REAL(wp), DIMENSION(jpij) :: zvti ! sum(v_i)
!
REAL(wp), DIMENSION(jpij,nlay_s) :: esrft ! snow energy of rafting ice
@@ -518,5 +515,7 @@
INTEGER , DIMENSION(jpij) :: itest_rdg, itest_rft ! test for conservation
!!-------------------------------------------------------------------
-
+ !
+ zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 ) ! total ice volume
+ !
! 1) Change in open water area due to closing and opening
!--------------------------------------------------------
@@ -535,6 +534,8 @@
IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! only if ice is ridging
- z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
-
+ IF( a_i_2d(ji,jl1) > epsi20 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
+ ELSE ; z1_ai(ji) = 0._wp
+ ENDIF
+
! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
@@ -549,10 +550,13 @@
! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
- vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg
+ 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?)
! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
virdg1 = v_i_2d (ji,jl1) * afrdg
- virdg2(ji) = v_i_2d (ji,jl1) * afrdg * ( 1. + rn_porordg )
+ virdg2(ji) = v_i_2d (ji,jl1) * afrdg + vsw
vsrdg(ji) = v_s_2d (ji,jl1) * afrdg
sirdg1 = sv_i_2d(ji,jl1) * afrdg
@@ -726,7 +730,8 @@
END DO ! jl1
!
+ ! roundoff errors
+ !----------------
! In case ridging/rafting lead to very small negative values (sometimes it happens)
- WHERE( a_i_2d(1:npti,:) < 0._wp ) a_i_2d(1:npti,:) = 0._wp
- WHERE( v_i_2d(1:npti,:) < 0._wp ) v_i_2d(1:npti,:) = 0._wp
+ 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, ze_s_2d, ze_i_2d )
!
END SUBROUTINE rdgrft_shift
@@ -854,5 +859,5 @@
CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) )
-
+ !
! !---------------------!
CASE( 2 ) !== from 1D to 2D ==!
@@ -945,5 +950,19 @@
CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
ENDIF
- ! ! allocate tke arrays
+ !
+ IF( .NOT. ln_icethd ) THEN
+ rn_porordg = 0._wp
+ rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp
+ rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp
+ IF( lwp ) THEN
+ WRITE(numout,*) ' ==> only ice dynamics is activated, thus some parameters must be changed'
+ WRITE(numout,*) ' rn_porordg = ', rn_porordg
+ WRITE(numout,*) ' rn_fsnwrdg = ', rn_fsnwrdg
+ WRITE(numout,*) ' rn_fpndrdg = ', rn_fpndrdg
+ WRITE(numout,*) ' rn_fsnwrft = ', rn_fsnwrft
+ WRITE(numout,*) ' rn_fpndrft = ', rn_fpndrft
+ ENDIF
+ ENDIF
+ ! ! allocate arrays
IF( ice_dyn_rdgrft_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
!
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rhg.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rhg.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rhg.F90 (revision 11083)
@@ -58,4 +58,6 @@
!!--------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ice time step
+ !
+ INTEGER :: jl ! dummy loop indices
!!--------------------------------------------------------------------
! controls
@@ -68,8 +70,15 @@
WRITE(numout,*)'~~~~~~~~~~~'
ENDIF
-
- ! --------
- ! Rheology
- ! --------
+ !
+ IF( ln_landfast_home ) THEN !-- Landfast ice parameterization
+ tau_icebfr(:,:) = 0._wp
+ DO jl = 1, jpl
+ WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
+ END DO
+ ENDIF
+ !
+ !--------------!
+ !== Rheology ==!
+ !--------------!
SELECT CASE( nice_rhg )
! !------------------------!
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rhg_evp.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rhg_evp.F90 (revision 11083)
@@ -112,5 +112,5 @@
REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i !
!!
- LOGICAL, PARAMETER :: ll_bdy_substep = .FALSE. ! temporary option to call bdy at each sub-time step (T)
+ LOGICAL, PARAMETER :: ll_bdy_substep = .TRUE. ! temporary option to call bdy at each sub-time step (T)
! or only at the main time step (F)
INTEGER :: ji, jj ! dummy loop indices
@@ -160,5 +160,6 @@
REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter
- REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity
+ REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small
+ REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small
!! --- diags
REAL(wp), DIMENSION(jpi,jpj) :: zswi
@@ -319,6 +320,8 @@
! switches
- zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
- zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
+ IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zswitchU(ji,jj) = 0._wp
+ ELSE ; zswitchU(ji,jj) = 1._wp ; ENDIF
+ IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; zswitchV(ji,jj) = 0._wp
+ ELSE ; zswitchV(ji,jj) = 1._wp ; ENDIF
END DO
@@ -519,5 +522,5 @@
& ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
- & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
+ & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
& ) * zmaskV(ji,jj)
ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
@@ -526,5 +529,5 @@
& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
- & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
+ & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
& ) * zmaskV(ji,jj)
ENDIF
@@ -567,5 +570,5 @@
& ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
- & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
+ & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
& ) * zmaskU(ji,jj)
ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
@@ -574,5 +577,5 @@
& ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
- & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
+ & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
& ) * zmaskU(ji,jj)
ENDIF
@@ -617,5 +620,5 @@
& ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
- & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
+ & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
& ) * zmaskU(ji,jj)
ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
@@ -624,5 +627,5 @@
& ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
- & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
+ & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
& ) * zmaskU(ji,jj)
ENDIF
@@ -665,5 +668,5 @@
& ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
- & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
+ & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
& ) * zmaskV(ji,jj)
ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
@@ -672,5 +675,5 @@
& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
- & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
+ & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
& ) * zmaskV(ji,jj)
ENDIF
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/iceitd.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/iceitd.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/iceitd.F90 (revision 11083)
@@ -21,4 +21,5 @@
USE ice1D ! sea-ice: thermodynamic variables
USE ice ! sea-ice: variables
+ USE icevar ! sea-ice: operations
USE icectl ! sea-ice: conservation tests
USE icetab ! sea-ice: convert 1D<=>2D
@@ -91,4 +92,6 @@
! 1) Identify grid cells with ice
!-----------------------------------------------------------------------------------------------
+ at_i(:,:) = SUM( a_i, dim=3 )
+ !
npti = 0 ; nptidx(:) = 0
DO jj = 1, jpj
@@ -249,5 +252,5 @@
! --- g(h) for each thickness category --- !
CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti) , a_i_1d(1:npti) , & ! in
- & g0 (1:npti,jl ), g1 (1:npti,jl), hL (1:npti,jl), hR (1:npti,jl) ) ! out
+ & g0 (1:npti,jl ), g1 (1:npti,jl), hL (1:npti,jl), hR (1:npti,jl) ) ! out
!
END DO
@@ -389,9 +392,11 @@
REAL(wp), DIMENSION(:,:), INTENT(in) :: pdvice ! ice volume transferred across boundary
!
- INTEGER :: ji, jj, jl, jk ! dummy loop indices
- INTEGER :: ii, ij, jl2, jl1 ! local integers
+ 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), DIMENSION(jpij,jpl) :: zaTsfn ! - -
+ REAL(wp), DIMENSION(jpij) :: zworka, zworkv ! 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
!!------------------------------------------------------------------
@@ -405,4 +410,14 @@
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), 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
+ ! to correct roundoff errors on a_i
+ CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d )
!----------------------------------------------------------------------------------------------
@@ -435,8 +450,4 @@
ELSE ; zworka(ji) = 0._wp
ENDIF
- !
- ! 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 )
- ! I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)
!
a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) ! Ice areas
@@ -476,8 +487,5 @@
!
DO jk = 1, nlay_s !--- Snow heat content
- !
DO ji = 1, npti
- ii = MOD( nptidx(ji) - 1, jpi ) + 1
- ij = ( nptidx(ji) - 1 ) / jpi + 1
!
jl1 = kdonor(ji,jl)
@@ -487,8 +495,7 @@
ELSE ; jl2 = jl
ENDIF
- !
- ztrans = e_s(ii,ij,jk,jl1) * zworkv(ji)
- e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
- e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
+ 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
@@ -497,6 +504,4 @@
DO jk = 1, nlay_i !--- Ice heat content
DO ji = 1, npti
- ii = MOD( nptidx(ji) - 1, jpi ) + 1
- ij = ( nptidx(ji) - 1 ) / jpi + 1
!
jl1 = kdonor(ji,jl)
@@ -506,8 +511,7 @@
ELSE ; jl2 = jl
ENDIF
- !
- ztrans = e_i(ii,ij,jk,jl1) * zworkv(ji)
- e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
- e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
+ 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
@@ -515,7 +519,21 @@
!
END DO ! boundaries, 1 to jpl-1
+
+ !-------------------
+ ! 3) roundoff errors
+ !-------------------
+ ! 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, ze_s_2d, ze_i_2d )
+
+ ! at_i must be <= rn_amax
+ zworka(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)
+ END DO
!-------------------------------------------------------------------------------
- ! 3) Update ice thickness and temperature
+ ! 4) Update ice thickness and temperature
!-------------------------------------------------------------------------------
WHERE( a_i_2d(1:npti,:) >= epsi20 )
@@ -536,4 +554,12 @@
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), 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
!
END SUBROUTINE itd_shiftice
@@ -558,4 +584,6 @@
!
IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution'
+ !
+ IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
!
jdonor(:,:) = 0
@@ -635,4 +663,6 @@
END DO
!
+ IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+ !
END SUBROUTINE ice_itd_reb
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icestp.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icestp.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icestp.F90 (revision 11083)
@@ -189,5 +189,5 @@
IF( ln_icethd ) CALL ice_thd( kt ) ! -- Ice thermodynamics
!
- IF( ln_icethd ) CALL ice_cor( kt , 2 ) ! -- Corrections
+ CALL ice_cor( kt , 2 ) ! -- Corrections
!
CALL ice_var_glo2eqv ! necessary calls (at least for coupling)
@@ -323,4 +323,7 @@
WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s
ENDIF
+ ! !--- change max ice concentration for roundoff errors
+ rn_amax_n = MIN( rn_amax_n, 1._wp - epsi10 )
+ rn_amax_s = MIN( rn_amax_s, 1._wp - epsi10 )
! !--- check consistency
IF ( jpl > 1 .AND. ln_virtual_itd ) THEN
@@ -431,7 +434,9 @@
t_si (:,:,:) = rt0 ! temp at the ice-snow interface
- tau_icebfr(:,:) = 0._wp ! landfast ice param only (clem: important to keep the init here)
- cnd_ice (:,:,:) = 0._wp ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T)
- qtr_ice_bot(:,:,:) = 0._wp ! initialization: part of solar radiation transmitted through the ice needed at least for outputs
+ tau_icebfr (:,:) = 0._wp ! landfast ice param only (clem: important to keep the init here)
+ cnd_ice (:,:,:) = 0._wp ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T)
+ qcn_ice (:,:,:) = 0._wp ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T)
+ qtr_ice_bot(:,:,:) = 0._wp ! initialization: part of solar radiation transmitted through the ice needed at least for outputs
+ qsb_ice_bot(:,:) = 0._wp ! (needed if ln_icethd=F)
!
! for control checks (ln_icediachk)
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd.F90 (revision 11083)
@@ -102,6 +102,4 @@
ENDIF
- CALL ice_var_glo2eqv
-
!---------------------------------------------!
! computation of friction velocity at T points
@@ -162,6 +160,7 @@
qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
- ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting
- IF( zqld > 0._wp ) THEN
+ ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting
+ ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget
+ IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN
fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
qlead(ji,jj) = 0._wp
@@ -178,7 +177,6 @@
! In case we bypass open-water ice formation
IF( .NOT. ln_icedO ) qlead(:,:) = 0._wp
- ! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere
+ ! In case we bypass growing/melting from top and bottom
IF( .NOT. ln_icedH ) THEN
- qt_atm_oi (:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
qsb_ice_bot(:,:) = 0._wp
fhld (:,:) = 0._wp
@@ -221,36 +219,32 @@
dh_i_sub (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp
dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp
- !
- IF( ln_icedH ) THEN ! --- growing/melting --- !
- CALL ice_thd_zdf ! Ice/Snow Temperature profile
- CALL ice_thd_dh ! Ice/Snow thickness
- CALL ice_thd_pnd ! Melt ponds formation
- CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping
+ !
+ CALL ice_thd_zdf ! --- Ice-Snow temperature --- !
+ !
+ IF( ln_icedH ) THEN ! --- Growing/Melting --- !
+ CALL ice_thd_dh ! Ice-Snow thickness
+ CALL ice_thd_pnd ! Melt ponds formation
+ CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping
ENDIF
- !
CALL ice_thd_sal( ln_icedS ) ! --- Ice salinity --- !
!
- CALL ice_thd_temp ! --- temperature update --- !
+ CALL ice_thd_temp ! --- Temperature update --- !
!
IF( ln_icedH .AND. ln_virtual_itd ) &
- & CALL ice_thd_mono ! --- extra lateral melting if virtual_itd --- !
- !
- IF( ln_icedA ) CALL ice_thd_da ! --- lateral melting --- !
+ & CALL ice_thd_mono ! --- Extra lateral melting if virtual_itd --- !
+ !
+ IF( ln_icedA ) CALL ice_thd_da ! --- Lateral melting --- !
!
CALL ice_thd_1d2d( jl, 2 ) ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
! ! --- & Move to 2D arrays --- !
- !
ENDIF
!
END DO
-
+ !
IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
- !
- CALL ice_var_zapsmall ! --- remove very small ice concentration (<1e-10) --- !
- ! ! & make sure at_i=SUM(a_i) & ato_i=1 where at_i=0
!
- IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- !
- !
- IF( ln_icedO ) CALL ice_thd_do ! --- frazil ice growing in leads --- !
+ IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- !
+ !
+ IF( ln_icedO ) CALL ice_thd_do ! --- Frazil ice growth in leads --- !
!
! controls
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd_do.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd_do.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd_do.F90 (revision 11083)
@@ -114,7 +114,5 @@
IF( ln_icediachk ) CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
- CALL ice_var_agg(1)
- CALL ice_var_glo2eqv
-
+ at_i(:,:) = SUM( a_i, dim=3 )
!------------------------------------------------------------------------------!
! 1) Collection thickness of ice formed in leads and polynyas
@@ -319,12 +317,12 @@
! --- lateral ice growth --- !
- ! If lateral ice growth gives an ice concentration gt 1, then
+ ! If lateral ice growth gives an ice concentration > amax, then
! we keep the excessive volume in memory and attribute it later to bottom accretion
DO ji = 1, npti
- IF ( za_newice(ji) > ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN
- zda_res(ji) = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) )
+ IF ( za_newice(ji) > MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error
+ zda_res(ji) = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) )
zdv_res(ji) = zda_res (ji) * zh_newice(ji)
- za_newice(ji) = za_newice(ji) - zda_res (ji)
- zv_newice(ji) = zv_newice(ji) - zdv_res (ji)
+ za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res (ji) )
+ zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res (ji) )
ELSE
zda_res(ji) = 0._wp
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd_zdf_bl99.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd_zdf_bl99.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd_zdf_bl99.F90 (revision 11083)
@@ -206,11 +206,11 @@
!
l_T_converged(:) = .FALSE.
- ! !============================!
! Convergence calculated until all sub-domain grid points have converged
! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation)
! but values are not taken into account (results independant of MPI partitioning)
!
+ ! !============================!
DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max ) ! Iterative procedure begins !
- ! !============================!
+ ! !============================!
iconv = iconv + 1
!
@@ -742,5 +742,5 @@
zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) )
! t_i
- DO jk = 0, nlay_i
+ DO jk = 1, nlay_i
ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0
t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
@@ -856,4 +856,8 @@
t_i_1d (1:npti,:) = ztiold (1:npti,:)
qcn_ice_1d(1:npti) = qcn_ice_top_1d(1:npti)
+
+ !!clem
+ ! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input
+ !clem
ENDIF
!
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icevar.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icevar.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icevar.F90 (revision 11083)
@@ -44,5 +44,6 @@
!! ice_var_salprof1d : salinity profile in the ice 1D
!! ice_var_zapsmall : remove very small area and volume
- !! ice_var_zapneg : remove negative ice fields (to debug the advection scheme UM3-5)
+ !! ice_var_zapneg : remove negative ice fields
+ !! ice_var_roundoff : remove negative values arising from roundoff erros
!! ice_var_itd : convert 1-cat to jpl-cat
!! ice_var_itd2 : convert N-cat to jpl-cat
@@ -71,4 +72,5 @@
PUBLIC ice_var_zapsmall
PUBLIC ice_var_zapneg
+ PUBLIC ice_var_roundoff
PUBLIC ice_var_itd
PUBLIC ice_var_itd2
@@ -229,5 +231,5 @@
IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
!
- ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]
+ ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]
ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C]
! Conversion q(S,T) -> T (second order equation)
@@ -236,5 +238,5 @@
t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts
!
- ELSE !--- no ice
+ ELSE !--- no ice
t_i(ji,jj,jk,jl) = rt0
ENDIF
@@ -537,5 +539,5 @@
- SUBROUTINE ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+ SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_zapneg ***
@@ -543,6 +545,5 @@
!! ** Purpose : Remove negative sea ice fields and correct fluxes
!!-------------------------------------------------------------------
- INTEGER :: ji, jj, jl, jk ! dummy loop indices
- !
+ REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
@@ -555,9 +556,16 @@
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
- !!-------------------------------------------------------------------
- !
+ !
+ INTEGER :: ji, jj, jl, jk ! dummy loop indices
+ REAL(wp) :: z1_dt
+ !!-------------------------------------------------------------------
+ !
+ z1_dt = 1._wp / pdt
!
DO jl = 1, jpl !== loop over the categories ==!
!
+ ! make sure a_i=0 where v_i<=0
+ WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp
+
!----------------------------------------
! zap ice energy and send it to the ocean
@@ -566,6 +574,6 @@
DO jj = 1 , jpj
DO ji = 1 , jpi
- IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
- hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 >0
+ IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+ hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
pe_i(ji,jj,jk,jl) = 0._wp
ENDIF
@@ -577,6 +585,6 @@
DO jj = 1 , jpj
DO ji = 1 , jpi
- IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
- hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
+ IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+ hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
pe_s(ji,jj,jk,jl) = 0._wp
ENDIF
@@ -590,14 +598,14 @@
DO jj = 1 , jpj
DO ji = 1 , jpi
- IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
- wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice
+ IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+ wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
pv_i (ji,jj,jl) = 0._wp
ENDIF
- IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
- wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice
+ IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+ wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
pv_s (ji,jj,jl) = 0._wp
ENDIF
- IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
- sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice
+ IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+ sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
psv_i (ji,jj,jl) = 0._wp
ENDIF
@@ -616,4 +624,36 @@
END SUBROUTINE ice_var_zapneg
+
+ SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
+ !!-------------------------------------------------------------------
+ !! *** ROUTINE ice_var_roundoff ***
+ !!
+ !! ** Purpose : Remove negative sea ice values arising from roundoff errors
+ !!-------------------------------------------------------------------
+ REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_i ! ice concentration
+ REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_i ! ice volume
+ REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_s ! snw volume
+ REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_i ! salt content
+ REAL(wp), DIMENSION(:,:) , INTENT(inout) :: poa_i ! age content
+ REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
+ REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume
+ REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content
+ REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content
+ !!-------------------------------------------------------------------
+ !
+ WHERE( pa_i (1:npti,:) < 0._wp .AND. pa_i (1:npti,:) > -epsi10 ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0
+ WHERE( pv_i (1:npti,:) < 0._wp .AND. pv_i (1:npti,:) > -epsi10 ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0
+ WHERE( pv_s (1:npti,:) < 0._wp .AND. pv_s (1:npti,:) > -epsi10 ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0
+ WHERE( psv_i(1:npti,:) < 0._wp .AND. psv_i(1:npti,:) > -epsi10 ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0
+ WHERE( poa_i(1:npti,:) < 0._wp .AND. poa_i(1:npti,:) > -epsi10 ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0
+ WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0
+ WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0
+ IF ( ln_pnd_H12 ) THEN
+ WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0
+ WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0
+ ENDIF
+ !
+ END SUBROUTINE ice_var_roundoff
+
SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
@@ -650,6 +690,6 @@
INTEGER :: idim, i_fill, jl0
REAL(wp) :: zarg, zV, zconv, zdh, zdv
- REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
- REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
+ REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
+ REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
INTEGER , DIMENSION(4) :: itest
!!-------------------------------------------------------------------
@@ -780,5 +820,5 @@
!!
!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
- !! by removing 25% ice area from jlmin and jlmax (resp.)
+ !! by removing 25% ice area from jlmin and jlmax (resp.)
!!
!! 3) Expand the filling to the empty cat between jlmin and jlmax
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icewri.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icewri.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icewri.F90 (revision 11083)
@@ -145,5 +145,5 @@
! record presence of fast ice
- WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk00(:,:) == 1._wp ) ; zfast(:,:) = 1._wp
+ WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp
ELSEWHERE ; zfast(:,:) = 0._wp
END WHERE
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/NST/agrif_top_update.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/NST/agrif_top_update.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/NST/agrif_top_update.F90 (revision 11083)
@@ -138,5 +138,5 @@
DO ji=i1,i2
IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
- trb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &
+ trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn) &
& + atfp * ( tabres_child(ji,jj,jk,jn) &
& - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdy_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdy_oce.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdy_oce.F90 (revision 11083)
@@ -85,5 +85,5 @@
!
INTEGER :: nb_bdy !: number of open boundary sets
- INTEGER :: nb_jpk_bdy !: number of levels in the bdy data (set < 0 if consistent with planned run)
+ INTEGER, DIMENSION(jp_bdy) :: nb_jpk_bdy !: number of levels in the bdy data (set < 0 if consistent with planned run)
INTEGER, DIMENSION(jp_bdy) :: nn_rimwidth !: boundary rim width for Flow Relaxation Scheme
INTEGER :: nn_volctl !: = 0 the total volume will have the variability of the surface Flux E-P
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydta.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydta.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydta.F90 (revision 11083)
@@ -243,5 +243,5 @@
IF( ln_full_vel_array(jbdy) ) THEN
CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), &
- & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy, &
+ & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy(jbdy), &
& fvl=ln_full_vel_array(jbdy) )
ELSE
@@ -313,5 +313,5 @@
jend = jstart + dta%nread(1) - 1
CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), &
- & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, &
+ & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy(jbdy), &
& fvl=ln_full_vel_array(jbdy) )
ENDIF
@@ -446,5 +446,5 @@
NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s
#endif
- NAMELIST/nambdy_dta/ ln_full_vel, nb_jpk_bdy
+ NAMELIST/nambdy_dta/ ln_full_vel
!!---------------------------------------------------------------------------
!
@@ -508,9 +508,9 @@
! Read namelists
! --------------
- REWIND(numnam_ref)
REWIND(numnam_cfg)
jfld = 0
DO jbdy = 1, nb_bdy
IF( nn_dta(jbdy) == 1 ) THEN
+ REWIND(numnam_ref)
READ ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp )
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydyn2d.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydyn2d.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydyn2d.F90 (revision 11083)
@@ -187,5 +187,5 @@
! Use characteristics method instead
zflag = ABS(flagu)
- zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij)
+ zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(ii+NINT(flagu),ij)
pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1)
END DO
@@ -205,5 +205,5 @@
! Use characteristics method instead
zflag = ABS(flagv)
- zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1)
+ zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ij+NINT(flagv))
pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
END DO
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyice.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyice.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyice.F90 (revision 11083)
@@ -57,6 +57,7 @@
INTEGER :: jbdy ! BDY set index
!!----------------------------------------------------------------------
- !
- IF( ln_timing ) CALL timing_start('bdy_ice_thd')
+ ! controls
+ IF( ln_timing ) CALL timing_start('bdy_ice_thd') ! timing
+ IF( ln_icediachk ) CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
!
CALL ice_var_glo2eqv
@@ -78,6 +79,8 @@
CALL ice_var_agg(1)
!
- IF( ln_icectl ) CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
- IF( ln_timing ) CALL timing_stop('bdy_ice_thd')
+ ! controls
+ IF( ln_icediachk ) CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+ IF( ln_icectl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints
+ IF( ln_timing ) CALL timing_stop ('bdy_ice_thd') ! timing
!
END SUBROUTINE bdy_ice
@@ -148,8 +151,8 @@
jpbound = 0 ; ib = ji ; jb = jj
!
- IF( u_ice(ji+1,jj ) < 0. .AND. umask(ji-1,jj ,1) == 0. ) jpbound = 1 ; ib = ji+1 ; jb = jj
- IF( u_ice(ji-1,jj ) > 0. .AND. umask(ji+1,jj ,1) == 0. ) jpbound = 1 ; ib = ji-1 ; jb = jj
- IF( v_ice(ji ,jj+1) < 0. .AND. vmask(ji ,jj-1,1) == 0. ) jpbound = 1 ; ib = ji ; jb = jj+1
- IF( v_ice(ji ,jj-1) > 0. .AND. vmask(ji ,jj+1,1) == 0. ) jpbound = 1 ; ib = ji ; jb = jj-1
+ IF( u_ice(ji ,jj ) < 0. .AND. umask(ji-1,jj ,1) == 0. ) jpbound = 1 ; ib = ji+1
+ IF( u_ice(ji-1,jj ) > 0. .AND. umask(ji ,jj ,1) == 0. ) jpbound = 1 ; ib = ji-1
+ IF( v_ice(ji ,jj ) < 0. .AND. vmask(ji ,jj-1,1) == 0. ) jpbound = 1 ; jb = jj+1
+ IF( v_ice(ji ,jj-1) > 0. .AND. vmask(ji ,jj ,1) == 0. ) jpbound = 1 ; jb = jj-1
!
IF( nn_ice_dta(jbdy) == 0 ) jpbound = 0 ; ib = ji ; jb = jj ! case ice boundaries = initial conditions
@@ -306,5 +309,5 @@
! one of the two zmsk is always 0 (because of zflag)
zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice
- zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice
+ zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) ) ! 0 if no ice
!
! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then do not change u_ice)
@@ -329,5 +332,5 @@
! one of the two zmsk is always 0 (because of zflag)
zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice
- zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice
+ zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) ) ! 0 if no ice
!
! v_ice = v_ice of the adjacent grid point except if this grid point is ice-free (then do not change v_ice)
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyini.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyini.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyini.F90 (revision 11083)
@@ -140,5 +140,5 @@
INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points
CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid
- INTEGER :: com_east, com_west, com_south, com_north ! Flags for boundaries sending
+ INTEGER :: com_east, com_west, com_south, com_north, jpk_max ! Flags for boundaries sending
INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b ! Flags for boundaries receiving
INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4) ! Arrays for neighbours coordinates
@@ -397,5 +397,5 @@
IF(lwp) WRITE(numout,*)
ENDIF
- IF( nb_jpk_bdy > 0 ) THEN
+ IF( nb_jpk_bdy(ib_bdy) > 0 ) THEN
IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***'
ELSE
@@ -516,25 +516,16 @@
ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), &
& nbrdta(jpbdta, jpbgrd, nb_bdy) )
-
- IF( nb_jpk_bdy>0 ) THEN
- ALLOCATE( dta_global(jpbdtau, 1, nb_jpk_bdy) )
- ALLOCATE( dta_global_z(jpbdtau, 1, nb_jpk_bdy) )
- ALLOCATE( dta_global_dz(jpbdtau, 1, nb_jpk_bdy) )
- ELSE
- ALLOCATE( dta_global(jpbdtau, 1, jpk) )
- ALLOCATE( dta_global_z(jpbdtau, 1, jpk) ) ! needed ?? TODO
- ALLOCATE( dta_global_dz(jpbdtau, 1, jpk) )! needed ?? TODO
- ENDIF
+
+ jpk_max = MAXVAL(nb_jpk_bdy)
+ jpk_max = MAX(jpk_max, jpk)
+
+ ALLOCATE( dta_global(jpbdtau, 1, jpk_max) )
+ ALLOCATE( dta_global_z(jpbdtau, 1, jpk_max) ) ! needed ?? TODO
+ ALLOCATE( dta_global_dz(jpbdtau, 1, jpk_max) )! needed ?? TODO
IF ( icount>0 ) THEN
- IF( nb_jpk_bdy>0 ) THEN
- ALLOCATE( dta_global2(jpbdtas, nrimmax, nb_jpk_bdy) )
- ALLOCATE( dta_global2_z(jpbdtas, nrimmax, nb_jpk_bdy) )
- ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, nb_jpk_bdy) )
- ELSE
- ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) )
- ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk) ) ! needed ?? TODO
- ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk) )! needed ?? TODO
- ENDIF
+ ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk_max) )
+ ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk_max) ) ! needed ?? TODO
+ ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk_max) )! needed ?? TODO
ENDIF
!
@@ -960,5 +951,5 @@
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
- if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
+ if( ii == (nlcit(nowe+1)-1) ) then
ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
@@ -974,5 +965,5 @@
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
- if((com_east_b .ne. 1) .and. (ii == 2)) then
+ if( ii == 2 ) then
ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
@@ -989,5 +980,5 @@
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
- if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
+ if( ii == (nlcit(nowe+1)-1) ) then
ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
@@ -1004,5 +995,5 @@
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
- if((com_east_b .ne. 1) .and. (ii == 2)) then
+ if( ii == 2 ) then
ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DIA/diacfl.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DIA/diacfl.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DIA/diacfl.F90 (revision 11083)
@@ -29,9 +29,4 @@
REAL(wp) :: rCu_max, rCv_max, rCw_max ! associated run max Courant number
-!!gm CAUTION: need to declare these arrays here, otherwise the calculation fails in multi-proc !
-!!gm I don't understand why.
- REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zCu_cfl, zCv_cfl, zCw_cfl ! workspace
-!!gm end
-
PUBLIC dia_cfl ! routine called by step.F90
PUBLIC dia_cfl_init ! routine called by nemogcm
@@ -55,8 +50,8 @@
INTEGER, INTENT(in) :: kt ! ocean time-step index
!
- INTEGER :: ji, jj, jk ! dummy loop indices
- REAL(wp) :: z2dt, zCu_max, zCv_max, zCw_max ! local scalars
- INTEGER , DIMENSION(3) :: iloc_u , iloc_v , iloc_w , iloc ! workspace
-!!gm this does not work REAL(wp), DIMENSION(jpi,jpj,jpk) :: zCu_cfl, zCv_cfl, zCw_cfl ! workspace
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: z2dt, zCu_max, zCv_max, zCw_max ! local scalars
+ INTEGER , DIMENSION(3) :: iloc_u , iloc_v , iloc_w , iloc ! workspace
+ REAL(wp), DIMENSION(jpi,jpj,jpk) :: zCu_cfl, zCv_cfl, zCw_cfl ! workspace
!!----------------------------------------------------------------------
!
@@ -71,5 +66,5 @@
DO jk = 1, jpk ! calculate Courant numbers
DO jj = 1, jpj
- DO ji = 1, fs_jpim1 ! vector opt.
+ DO ji = 1, jpi
zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u (ji,jj) ! for i-direction
zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v (ji,jj) ! for j-direction
@@ -111,5 +106,5 @@
! ! write out to file
IF( lwp ) THEN
- WRITE(numcfl,FMT='(2x,i4,5x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
+ WRITE(numcfl,FMT='(2x,i6,3x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
WRITE(numcfl,FMT='(11x, a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') 'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3)
WRITE(numcfl,FMT='(11x, a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') 'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3)
@@ -172,8 +167,4 @@
rCw_max = 0._wp
!
-!!gm required to work
- ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) )
-!!gm end
- !
END SUBROUTINE dia_cfl_init
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/dynkeg.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/dynkeg.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/dynkeg.F90 (revision 11083)
@@ -74,10 +74,10 @@
INTEGER, INTENT( in ) :: kscheme ! =0/1 type of KEG scheme
!
- INTEGER :: ji, jj, jk, jb ! dummy loop indices
- INTEGER :: ii, ifu, ib_bdy ! local integers
- INTEGER :: ij, ifv, igrd ! - -
- REAL(wp) :: zu, zv ! local scalars
+ INTEGER :: ji, jj, jk, jb ! dummy loop indices
+ INTEGER :: ifu, ifv, igrd, ib_bdy ! local integers
+ REAL(wp) :: zu, zv ! local scalars
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhke
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv
+ REAL(wp) :: zweightu, zweightv
!!----------------------------------------------------------------------
!
@@ -97,32 +97,4 @@
zhke(:,:,jpk) = 0._wp
-
- IF (ln_bdy) THEN
- ! Maria Luneva & Fred Wobus: July-2016
- ! compensate for lack of turbulent kinetic energy on liquid bdy points
- DO ib_bdy = 1, nb_bdy
- IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
- igrd = 2 ! Copying normal velocity into points outside bdy
- DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
- DO jk = 1, jpkm1
- ii = idx_bdy(ib_bdy)%nbi(jb,igrd)
- ij = idx_bdy(ib_bdy)%nbj(jb,igrd)
- ifu = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
- un(ii-ifu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)
- END DO
- END DO
- !
- igrd = 3 ! Copying normal velocity into points outside bdy
- DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
- DO jk = 1, jpkm1
- ii = idx_bdy(ib_bdy)%nbi(jb,igrd)
- ij = idx_bdy(ib_bdy)%nbj(jb,igrd)
- ifv = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
- vn(ii,ij-ifv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)
- END DO
- END DO
- ENDIF
- ENDDO
- ENDIF
SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==!
@@ -140,4 +112,28 @@
END DO
END DO
+ !
+ IF (ln_bdy) THEN
+ ! Maria Luneva & Fred Wobus: July-2016
+ ! compensate for lack of turbulent kinetic energy on liquid bdy points
+ DO ib_bdy = 1, nb_bdy
+ IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
+ igrd = 1 ! compensating null velocity on the bdy
+ DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+ ji = idx_bdy(ib_bdy)%nbi(jb,igrd) ! maximum extent : from 2 to jpi-1
+ jj = idx_bdy(ib_bdy)%nbj(jb,igrd) ! maximum extent : from 2 to jpj-1
+ DO jk = 1, jpkm1
+ zhke(ji,jj,jk) = 0._wp
+ zweightu = umask(ji-1,jj ,jk) + umask(ji,jj,jk)
+ zweightv = vmask(ji ,jj-1,jk) + vmask(ji,jj,jk)
+ zu = un(ji-1,jj ,jk) * un(ji-1,jj ,jk) + un(ji ,jj ,jk) * un(ji ,jj ,jk)
+ zv = vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) + vn(ji ,jj ,jk) * vn(ji ,jj ,jk)
+ IF( zweightu > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zu / (2._wp * zweightu)
+ IF( zweightv > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zv / (2._wp * zweightv)
+ END DO
+ END DO
+ END IF
+ CALL lbc_bdy_lnk( 'dynkeg', zhke, 'T', 1., ib_bdy ) ! send 2 and recv jpi, jpj used in the computation of the speed tendencies
+ END DO
+ END IF
!
CASE ( nkeg_HW ) !-- Hollingsworth scheme --!
@@ -158,14 +154,37 @@
END DO
END DO
+ IF (ln_bdy) THEN
+ ! Maria Luneva & Fred Wobus: July-2016
+ ! compensate for lack of turbulent kinetic energy on liquid bdy points
+ DO ib_bdy = 1, nb_bdy
+ IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
+ igrd = 1 ! compensation null velocity on land at the bdy
+ DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+ ji = idx_bdy(ib_bdy)%nbi(jb,igrd) ! maximum extent : from 2 to jpi-1
+ jj = idx_bdy(ib_bdy)%nbj(jb,igrd) ! maximum extent : from 2 to jpj-1
+ DO jk = 1, jpkm1
+ zhke(ji,jj,jk) = 0._wp
+ zweightu = 8._wp * ( umask(ji-1,jj ,jk) + umask(ji ,jj ,jk) ) &
+ & + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji ,jj-1,jk) + umask(ji ,jj+1,jk) )
+ zweightv = 8._wp * ( vmask(ji ,jj-1,jk) + vmask(ji ,jj-1,jk) ) &
+ & + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj ,jk) + vmask(ji+1,jj ,jk) )
+ zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) &
+ & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) &
+ & + ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) &
+ & + ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) * ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) )
+ zv = 8._wp * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) &
+ & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) &
+ & + ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) &
+ & + ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) * ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) )
+ IF( zweightu > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zu / ( 2._wp * zweightu )
+ IF( zweightv > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zv / ( 2._wp * zweightv )
+ END DO
+ END DO
+ END IF
+ END DO
+ END IF
CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
!
- END SELECT
-
- IF (ln_bdy) THEN
- ! restore velocity masks at points outside boundary
- un(:,:,:) = un(:,:,:) * umask(:,:,:)
- vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)
- ENDIF
-
+ END SELECT
!
DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==!
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/sshwzv.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/sshwzv.F90 (revision 11083)
@@ -297,7 +297,6 @@
IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
IF(lwp) WRITE(numout,*) '~~~~~ '
- !
- Cu_adv(:,:,jpk) = 0._wp ! bottom value : Cu_adv=0 (set once for all)
- ENDIF
+ ENDIF
+ !
!
DO jk = 1, jpkm1 ! calculate Courant numbers
@@ -305,27 +304,27 @@
DO ji = 2, fs_jpim1 ! vector opt.
z1_e3w = 1._wp / e3w_n(ji,jj,jk)
- Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) &
- & + ( MAX( e2u(ji ,jj)*e3uw_n(ji ,jj,jk)*un(ji ,jj,jk), 0._wp ) - &
- & MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) ) &
- & * r1_e1e2t(ji,jj) &
- & + ( MAX( e1v(ji,jj )*e3vw_n(ji,jj ,jk)*vn(ji,jj ,jk), 0._wp ) - &
- & MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) ) &
- & * r1_e1e2t(ji,jj) &
- & ) * z1_e3w
+ Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) & ! 2*rdt and not r2dt (for restartability)
+ & + ( MAX( e2u(ji ,jj)*e3uw_n(ji ,jj,jk)*un(ji ,jj,jk), 0._wp ) - &
+ & MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) ) &
+ & * r1_e1e2t(ji,jj) &
+ & + ( MAX( e1v(ji,jj )*e3vw_n(ji,jj ,jk)*vn(ji,jj ,jk), 0._wp ) - &
+ & MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) ) &
+ & * r1_e1e2t(ji,jj) &
+ & ) * z1_e3w
END DO
END DO
END DO
+ CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. )
!
CALL iom_put("Courant",Cu_adv)
!
- wi(:,:,:) = 0._wp ! Includes top and bottom values set to zero
IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere
DO jk = 1, jpkm1 ! or scan Courant criterion and partition
- DO jj = 2, jpjm1 ! w where necessary
- DO ji = 2, fs_jpim1 ! vector opt.
+ DO jj = 1, jpj ! w where necessary
+ DO ji = 1, jpi
!
zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
!
- IF( zCu < Cu_min ) THEN !<-- Fully explicit
+ IF( zCu <= Cu_min ) THEN !<-- Fully explicit
zcff = 0._wp
ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit
@@ -346,5 +345,6 @@
ELSE
! Fully explicit everywhere
- Cu_adv = 0.0_wp ! Reuse array to output coefficient
+ Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient
+ wi (:,:,:) = 0._wp
ENDIF
CALL iom_put("wimp",wi)
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/LBC/lib_mpp.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/LBC/lib_mpp.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/LBC/lib_mpp.F90 (revision 11083)
@@ -1480,6 +1480,7 @@
LOGICAL , OPTIONAL, INTENT(in ) :: ld_lbc, ld_glb, ld_dlg
!!
+ CHARACTER(len=128) :: ccountname ! name of a subroutine to count communications
LOGICAL :: ll_lbc, ll_glb, ll_dlg
- INTEGER :: ji, jj, jk, jh, jf ! dummy loop indices
+ INTEGER :: ji, jj, jk, jh, jf, jcount ! dummy loop indices
!!----------------------------------------------------------------------
!
@@ -1538,13 +1539,21 @@
WRITE(numcom,*) ' '
WRITE(numcom,*) ' lbc_lnk called'
- jj = 1
- DO ji = 2, n_sequence_lbc
- IF( crname_lbc(ji-1) /= crname_lbc(ji) ) THEN
- WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_lbc(ji-1))
- jj = 0
+ DO ji = 1, n_sequence_lbc - 1
+ IF ( crname_lbc(ji) /= 'already counted' ) THEN
+ ccountname = crname_lbc(ji)
+ crname_lbc(ji) = 'already counted'
+ jcount = 1
+ DO jj = ji + 1, n_sequence_lbc
+ IF ( ccountname == crname_lbc(jj) ) THEN
+ jcount = jcount + 1
+ crname_lbc(jj) = 'already counted'
+ END IF
+ END DO
+ WRITE(numcom,'(A, I4, A, A)') ' - ', jcount,' times by subroutine ', TRIM(ccountname)
END IF
- jj = jj + 1
END DO
- WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_lbc(n_sequence_lbc))
+ IF ( crname_lbc(n_sequence_lbc) /= 'already counted' ) THEN
+ WRITE(numcom,'(A, I4, A, A)') ' - ', 1,' times by subroutine ', TRIM(crname_lbc(ncom_rec_max))
+ END IF
WRITE(numcom,*) ' '
IF ( n_sequence_glb > 0 ) THEN
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/ZDF/zdfphy.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/ZDF/zdfphy.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/ZDF/zdfphy.F90 (revision 11083)
@@ -132,6 +132,7 @@
IF( ln_zad_Aimp ) THEN
IF( zdf_phy_alloc() /= 0 ) &
- & CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' )
- wi(:,:,:) = 0._wp
+ & CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' )
+ Cu_adv(:,:,:) = 0._wp
+ wi (:,:,:) = 0._wp
ENDIF
! !== Background eddy viscosity and diffusivity ==!
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/TOP/PISCES/P4Z/p4zpoc.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/TOP/PISCES/P4Z/p4zpoc.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/TOP/PISCES/P4Z/p4zpoc.F90 (revision 11083)
@@ -166,5 +166,5 @@
& + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1) / tgfunc(ji,jj,jk-1) * ( 1. &
& - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn) * zsizek ) + prodgoc(ji,jj,jk) &
- & / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn)
+ & / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn) * alphan(jn)
alphat = alphat + alphag(ji,jj,jk,jn)
remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90 (revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90 (revision 11083)
@@ -96,8 +96,8 @@
! ! ==>>> set by hand non-zero value on first/last columns & rows
DO ji = mi0(1), mi1(1) ! first row of global domain only
- zhu(ji,2) = zht(1,2)
- END DO
- DO ji = mi0(jpi), mi1(jpi) ! last row of global domain only
- zhu(ji,2) = zht(jpi,2)
+ zhu(ji,2) = zht(ji,2)
+ END DO
+ DO ji = mi0(jpiglo), mi1(jpiglo) ! last row of global domain only
+ zhu(ji,2) = zht(ji,2)
END DO
zhu(:,1) = zhu(:,2)