Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 2154 additions and 284 deletions
src/OCE/TRA/traqsr.F90
View file @ d1c6ae27
......@@ -12,10 +12,16 @@ MODULE traqsr
!! 3.6 ! 2012-05 (C. Rousset) store attenuation coef for use in ice model
!! 3.6 ! 2015-12 (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll
!! 3.7 ! 2015-11 (G. Madec, A. Coward) remove optimisation for fix volume
!! 4.0 ! 2020-11 (A. Coward) optimisation
!! 4.5 ! 2021-03 (G. Madec) further optimisation + adaptation for RK3
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! tra_qsr : temperature trend due to the penetration of solar radiation
!! qsr_RGBc : IR + RGB light penetration with Chlorophyll data case
!! qsr_RGB : IR + RGB light penetration with constant Chlorophyll case
!! qsr_2BD : 2 bands (InfraRed + Visible light) case
!! qsr_ext_lev : level of extinction for each bands
!! tra_qsr_init : initialization of the qsr penetration
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and active tracers
......@@ -52,16 +58,25 @@ MODULE traqsr
REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands)
REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands)
!
INTEGER , PUBLIC :: nksr !: levels below which the light cannot penetrate (depth larger than 391 m)
INTEGER, PARAMETER :: np_RGB = 1 ! R-G-B light penetration with constant Chlorophyll
INTEGER, PARAMETER :: np_RGBc = 2 ! R-G-B light penetration with Chlorophyll data
INTEGER, PARAMETER :: np_2BD = 3 ! 2 bands light penetration
INTEGER, PARAMETER :: np_BIO = 4 ! bio-model light penetration
!
INTEGER :: nqsr ! user choice of the type of light penetration
REAL(wp) :: xsi0r ! inverse of rn_si0
REAL(wp) :: xsi1r ! inverse of rn_si1
INTEGER :: nqsr ! user choice of the type of light penetration
INTEGER :: nc_rgb ! RGB with cst Chlorophyll: index associated with the chosen Chl value
!
! ! extinction level
INTEGER :: nk0 !: IR (depth larger ~12 m)
INTEGER :: nkV !: Visible light (depth larger than ~840 m)
INTEGER :: nkR, nkG, nkB !: RGB (depth larger than ~100 m, ~470 m, ~1700 m, resp.)
!
INTEGER, PUBLIC :: nksr !: =nkV, i.e. maximum level of light extinction (used in traatf(_qco).F90)
!
! ! inverse of attenuation length
REAL(wp) :: r1_si0 ! all schemes : infrared = 1/rn_si0
REAL(wp) :: r1_si1 ! 2 band : mean RGB = 1/rn_si1
REAL(wp) :: r1_LR, r1_LG, r1_LB ! RGB with constant Chl (np_RGB)
!
REAL(wp) , PUBLIC, DIMENSION(3,61) :: rkrgb ! tabulated attenuation coefficients for RGB absorption
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read)
......@@ -71,7 +86,7 @@ MODULE traqsr
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: traqsr.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: traqsr.F90 15157 2021-07-29 08:28:32Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -84,38 +99,20 @@ CONTAINS
!! penetration and add it to the general temperature trend.
!!
!! ** Method : The profile of the solar radiation within the ocean is defined
!! through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
!! Considering the 2 wavebands case:
!! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
!! The temperature trend associated with the solar radiation penetration
!! is given by : zta = 1/e3t dk[ I ] / (rho0*Cp)
!! At the bottom, boudary condition for the radiation is no flux :
!! all heat which has not been absorbed in the above levels is put
!! in the last ocean level.
!! through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) or computed by
!! the biogeochemical model
!! The computation is only done down to the level where
!! I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) .
!! I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) .
!!
!! ** Action : - update ta with the penetrative solar radiation trend
!! ** Action : - update ts(jp_tem) with the penetrative solar radiation trend
!! - send trend for further diagnostics (l_trdtra=T)
!!
!! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
!! Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
!! Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
!!----------------------------------------------------------------------
INTEGER, INTENT(in ) :: kt ! ocean time-step
INTEGER, INTENT(in ) :: Kmm, Krhs ! time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
INTEGER, INTENT(in ) :: kt, Kmm, Krhs ! ocean time-step and time-level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: irgb ! local integers
REAL(wp) :: zchl, zcoef, z1_2 ! local scalars
REAL(wp) :: zc0 , zc1 , zc2 , zc3 ! - -
REAL(wp) :: zzc0, zzc1, zzc2, zzc3 ! - -
REAL(wp) :: zz0 , zz1 , ze3t, zlui ! - -
REAL(wp) :: zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze
REAL(wp) :: zlogc, zlogze, zlogCtot, zlogCze
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ze0, ze1, ze2, ze3
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot, ztmp3d
REAL(wp) :: z1_2, ze3t ! local scalars
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('tra_qsr')
......@@ -128,10 +125,13 @@ CONTAINS
ENDIF
ENDIF
!
IF( l_trdtra ) THEN ! trends diagnostic: save the input temperature trend
IF( l_trdtra ) THEN ! trends diagnostic: save the input temperature trend
ALLOCATE( ztrdt(jpi,jpj,jpk) )
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
ENDIF
!
#if ! defined key_RK3
! ! MLF only : heat content trend due to Qsr flux (qsr_hc)
!
! !-----------------------------------!
! ! before qsr induced heat content !
......@@ -145,140 +145,62 @@ CONTAINS
ENDIF
ELSE ! No restart or Euler forward at 1st time step
z1_2 = 1._wp
DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
DO_3D_OVR( 0, 0, 0, 0, 1, jpk )
qsr_hc_b(ji,jj,jk) = 0._wp
END_3D
ENDIF
ELSE !== Swap of qsr heat content ==!
z1_2 = 0.5_wp
DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
DO_3D_OVR( 0, 0, 0, 0, 1, jpk )
qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk)
END_3D
ENDIF
!
! !--------------------------------!
SELECT CASE( nqsr ) ! now qsr induced heat content !
! !--------------------------------!
!
CASE( np_BIO ) !== bio-model fluxes ==!
!
DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr )
qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) )
END_3D
!
CASE( np_RGB , np_RGBc ) !== R-G-B fluxes ==!
!
ALLOCATE( ze0 (A2D(nn_hls)) , ze1 (A2D(nn_hls)) , &
& ze2 (A2D(nn_hls)) , ze3 (A2D(nn_hls)) , &
& ztmp3d(A2D(nn_hls),nksr + 1) )
!
IF( nqsr == np_RGBc ) THEN !* Variable Chlorophyll
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only for the full domain
IF( ln_tile ) CALL dom_tile_stop( ldhold=.TRUE. ) ! Use full domain
CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step
IF( ln_tile ) CALL dom_tile_start( ldhold=.TRUE. ) ! Revert to tile domain
ENDIF
!
! Separation in R-G-B depending on the surface Chl
! perform and store as many of the 2D calculations as possible
! before the 3D loop (use the temporary 2D arrays to replace the
! most expensive calculations)
!
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
! zlogc = log(zchl)
zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) )
! zc1 : log(zCze) = log (1.12 * zchl**0.803)
zc1 = 0.113328685307 + 0.803 * zlogc
! zc2 : log(zCtot) = log(40.6 * zchl**0.459)
zc2 = 3.703768066608 + 0.459 * zlogc
! zc3 : log(zze) = log(568.2 * zCtot**(-0.746))
zc3 = 6.34247346942 - 0.746 * zc2
! IF( log(zze) > log(102.) ) log(zze) = log(200.0 * zCtot**(-0.293))
IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2
!
ze0(ji,jj) = zlogc ! ze0 = log(zchl)
ze1(ji,jj) = EXP( zc1 ) ! ze1 = zCze
ze2(ji,jj) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) ) ! ze2 = 1/zdelpsi
ze3(ji,jj) = EXP( - zc3 ) ! ze3 = 1/zze
END_2D
#endif
!
DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr + 1 )
! zchl = ALOG( ze0(ji,jj) )
zlogc = ze0(ji,jj)
!
zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
!
zCze = ze1(ji,jj)
zrdpsi = ze2(ji,jj) ! 1/zdelpsi
zpsi = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze
!
! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) )
zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) )
! Convert chlorophyll value to attenuation coefficient look-up table index
ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15
END_3D
ELSE !* constant chlorophyll
zchl = 0.05
! NB. make sure constant value is such that:
zchl = MIN( 10. , MAX( 0.03, zchl ) )
! Convert chlorophyll value to attenuation coefficient look-up table index
zlui = 41 + 20.*LOG10(zchl) + 1.e-15
DO jk = 1, nksr + 1
ztmp3d(:,:,jk) = zlui
END DO
ENDIF
!
zcoef = ( 1. - rn_abs ) / 3._wp !* surface equi-partition in R-G-B
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
ze0(ji,jj) = rn_abs * qsr(ji,jj)
ze1(ji,jj) = zcoef * qsr(ji,jj)
ze2(ji,jj) = zcoef * qsr(ji,jj)
ze3(ji,jj) = zcoef * qsr(ji,jj)
! store the surface SW radiation; re-use the surface ztmp3d array
! since the surface attenuation coefficient is not used
ztmp3d(ji,jj,1) = qsr(ji,jj)
END_2D
!
! !* interior equi-partition in R-G-B depending on vertical profile of Chl
DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr + 1 )
ze3t = e3t(ji,jj,jk-1,Kmm)
irgb = NINT( ztmp3d(ji,jj,jk) )
zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r )
zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) )
zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) )
zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) )
ze0(ji,jj) = zc0
ze1(ji,jj) = zc1
ze2(ji,jj) = zc2
ze3(ji,jj) = zc3
ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
END_3D
!
DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr ) !* now qsr induced heat content
qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
END_3D
! !----------------------------!
SELECT CASE( nqsr ) ! qsr induced heat content !
! !----------------------------!
!
CASE( np_RGBc ) ; CALL qsr_RGBc( kt, Kmm, pts, Krhs ) !== R-G-B fluxes using chlorophyll data ==! with Morel &Berthon (1989) vertical profile
!
DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d )
CASE( np_RGB ) ; CALL qsr_RGB ( kt, Kmm, pts, Krhs ) !== R-G-B fluxes with constant chlorophyll ==!
!
CASE( np_2BD ) !== 2-bands fluxes ==!
CASE( np_2BD ) ; CALL qsr_2BD ( Kmm, pts, Krhs ) !== 2-bands fluxes ==!
!
zz0 = rn_abs * r1_rho0_rcp ! surface equi-partition in 2-bands
zz1 = ( 1. - rn_abs ) * r1_rho0_rcp
DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr ) !* now qsr induced heat content
zc0 = zz0 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi1r )
zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) )
CASE( np_BIO ) !== bio-model fluxes ==!
DO_3D( 0, 0, 0, 0, 1, nkV )
#if defined key_RK3
! !- RK3 : temperature trend at jk t-level
ze3t = e3t(ji,jj,jk,Kmm)
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / ze3t
#else
! !- MLF : heat content trend due to Qsr flux (qsr_hc)
qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) )
#endif
END_3D
! !- sea-ice : store the 1st level attenuation coefficient
WHERE( etot3(A2D(0),1) /= 0._wp ) ; fraqsr_1lev(A2D(0)) = 1._wp - etot3(A2D(0),2) / etot3(A2D(0),1)
ELSEWHERE ; fraqsr_1lev(A2D(0)) = 1._wp
END WHERE
!
END SELECT
!
! !-----------------------------!
! ! update to the temp. trend !
! !-----------------------------!
#if defined key_RK3
! ! RK3 : diagnostics/output
IF( l_trdtra .OR. iom_use('qsr3d') ) THEN ! qsr diagnostics
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
! ! qsr tracers trends saved for diagnostics
IF( l_trdtra ) CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
IF( iom_use('qsr3d') ) THEN ! qsr distribution
DO jk = nkV, 1, -1
ztrdt(:,:,jk) = ztrdt(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp
END DO
CALL iom_put( 'qsr3d', ztrdt ) ! 3D distribution of shortwave Radiation
ENDIF
DEALLOCATE( ztrdt )
ENDIF
#else
! ! MLF : add the temperature trend
DO_3D( 0, 0, 0, 0, 1, nksr )
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) &
& + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) &
......@@ -286,7 +208,7 @@ CONTAINS
END_3D
!
! sea-ice: store the 1st ocean level attenuation coefficient
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 0, 0, 0 )
IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
ELSE ; fraqsr_1lev(ji,jj) = 1._wp
ENDIF
......@@ -301,6 +223,13 @@ CONTAINS
CALL iom_put( 'qsr3d', zetot ) ! 3D distribution of shortwave Radiation
DEALLOCATE( zetot )
ENDIF
!
IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
DEALLOCATE( ztrdt )
ENDIF
#endif
!
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile
IF( lrst_oce ) THEN ! write in the ocean restart file
......@@ -309,11 +238,6 @@ CONTAINS
ENDIF
ENDIF
!
IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
DEALLOCATE( ztrdt )
ENDIF
! ! print mean trends (used for debugging)
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr - Ta: ', mask1=tmask, clinfo3='tra-ta' )
!
......@@ -322,6 +246,518 @@ CONTAINS
END SUBROUTINE tra_qsr
SUBROUTINE qsr_RGBc( kt, Kmm, pts, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE qsr_RGBc ***
!!
!! ** Purpose : Red-Green-Blue solar radiation using chlorophyll data
!!
!! ** Method : The profile of the solar radiation within the ocean is defined
!! through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
!! Considering the 2 wavebands case:
!! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
!! The temperature trend associated with the solar radiation penetration
!! is given by : zta = 1/e3t dk[ I ] / (rho0*Cp)
!! At the bottom, boudary condition for the radiation is no flux :
!! all heat which has not been absorbed in the above levels is put
!! in the last ocean level.
!! The computation is only done down to the level where
!! I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) .
!!
!! ** Action : - update ta with the penetrative solar radiation trend
!! - send trend for further diagnostics (l_trdtra=T)
!!
!! Reference : Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
!! Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
!!----------------------------------------------------------------------
INTEGER, INTENT(in ) :: kt, Kmm, Krhs ! ocean time-step and time-level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
!!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: irgb ! local integer
REAL(wp) :: zc1 , zc2 , zc3, zchl ! local scalars
REAL(wp) :: zze0, zzeR, zzeG, zzeB, zzeT ! - -
REAL(wp) :: zz0 , zz1 , ze3t ! - -
REAL(wp) :: zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze ! - -
REAL(wp) :: zlogc, zlogze, zlogCtot, zlogCze ! - -
REAL(wp), DIMENSION(A2D(0) ) :: ze0, zeR, zeG, zeB, zeT
REAL(wp), DIMENSION(A2D(0),0:3) :: zc
!!----------------------------------------------------------------------
!
!
! !===========================================!
! !== R-G-B fluxes using chlorophyll data ==! with Morel &Berthon (1989) vertical profile
! !===================================****====!
!
! != Chlorophyll data =!
!
IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only for the full domain
IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) ! Use full domain
CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step
IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 ) ! Revert to tile domain
ENDIF
!
DO_2D( 0, 0, 0, 0 ) ! pre-calculated expensive coefficient
zlogc = LOG( MAX( 0.03_wp, MIN( sf_chl(1)%fnow(ji,jj,1) ,10._wp ) ) ) ! zlogc = log(zchl) with 0.03 <= Chl >= 10.
zc1 = 0.113328685307 + 0.803 * zlogc ! zc1 : log(zCze) = log (1.12 * zchl**0.803)
zc2 = 3.703768066608 + 0.459 * zlogc ! zc2 : log(zCtot) = log(40.6 * zchl**0.459)
zc3 = 6.34247346942 - 0.746 * zc2 ! zc3 : log(zze) = log(568.2 * zCtot**(-0.746))
IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2 ! IF(log(zze)>log(102)) log(zze) = log(200*zCtot**(-0.293))
!
zc(ji,jj,0) = zlogc ! ze(0) = log(zchl)
zc(ji,jj,1) = EXP( zc1 ) ! ze(1) = zCze
zc(ji,jj,2) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) ) ! ze(2) = 1/zdelpsi
zc(ji,jj,3) = EXP( - zc3 ) ! ze(3) = 1/zze
END_2D
!
! != surface light =!
!
zz0 = rn_abs ! Infrared absorption
zz1 = ( 1._wp - rn_abs ) / 3._wp ! R-G-B equi-partition
!
DO_2D( 0, 0, 0, 0 ) ! surface light
ze0(ji,jj) = zz0 * qsr(ji,jj) ; zeR(ji,jj) = zz1 * qsr(ji,jj) ! IR ; Red
zeG(ji,jj) = zz1 * qsr(ji,jj) ; zeB(ji,jj) = zz1 * qsr(ji,jj) ! Green ; Blue
zeT(ji,jj) = qsr(ji,jj) ! Total
END_2D
!
! != interior light =!
!
DO jk = 1, nk0 !* near surface layers *! (< ~12 meters : IR + RGB )
DO_2D( 0, 0, 0, 0 )
! !- inverse of RGB attenuation lengths
zlogc = zc(ji,jj,0)
zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
zCze = zc(ji,jj,1)
zrdpsi = zc(ji,jj,2) ! 1/zdelpsi
!!st05 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze
zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm) ! gdepw/zze
! ! make sure zchl value is such that: 0.03 < zchl < 10.
zchl = MAX( 0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp ) )
! ! Convert chlorophyll value to attenuation coefficient
irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) ! look-up table index
! Red ! Green ! Blue
r1_LR = rkrgb(3,irgb) ; r1_LG = rkrgb(2,irgb) ; r1_LB = rkrgb(1,irgb)
!
! !- fluxes at jk+1 w-level
ze3t = e3t(ji,jj,jk,Kmm)
zze0 = ze0(ji,jj) * EXP( - ze3t*r1_si0 ) ; zzeR = zeR(ji,jj) * EXP( - ze3t*r1_LR ) ! IR ; Red at jk+1 w-level
zzeG = zeG(ji,jj) * EXP( - ze3t*r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t*r1_LB ) ! Green ; Blue - -
zzeT = ( zze0 + zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1) ! Total - -
!!st01 zzeT = ( zze0 + zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - -
!
#if defined key_RK3
! !- RK3 : temperature trend at jk t-level
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
#else
! !- MLF : heat content trend due to Qsr flux (qsr_hc)
qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
#endif
ze0(ji,jj) = zze0 ; zeR(ji,jj) = zzeR ! IR ; Red store at jk+1 w-level
zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue - - -
zeT(ji,jj) = zzeT ! total - - -
END_2D
!
END DO
!
DO jk = nk0+1, nkR !* down to Red extinction *! (< ~71 meters : RGB , IR removed from calculation)
DO_2D( 0, 0, 0, 0 )
! !- inverse of RGB attenuation lengths
zlogc = zc(ji,jj,0)
zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
zCze = zc(ji,jj,1)
zrdpsi = zc(ji,jj,2) ! 1/zdelpsi
zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm) ! gdepw/zze
!!st05 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze
! ! make sure zchl value is such that: 0.03 < zchl < 10.
zchl = MAX( 0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp ) )
! ! Convert chlorophyll value to attenuation coefficient
irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) ! look-up table index
! Red ! Green ! Blue
r1_LR = rkrgb(3,irgb) ; r1_LG = rkrgb(2,irgb) ; r1_LB = rkrgb(1,irgb)
!
! !- fluxes at jk+1 w-level
ze3t = e3t(ji,jj,jk,Kmm)
zzeR = zeR(ji,jj) * EXP( - ze3t*r1_LR ) ! Red at jk+1 w-level
zzeG = zeG(ji,jj) * EXP( - ze3t*r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t*r1_LB ) ! Green ; Blue - -
zzeT = ( zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - -
!
#if defined key_RK3
! !- RK3 : temperature trend at jk t-level
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
#else
! !- MLF : heat content trend due to Qsr flux (qsr_hc)
qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
#endif
zeR(ji,jj) = zzeR ! Red store at jk+1 w-level
zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue - - -
zeT(ji,jj) = zzeT ! total - - -
END_2D
END DO
!
DO jk = nkR+1, nkG !* down to Green extinction *! (< ~350 m : GB , IR+R removed from calculation)
DO_2D( 0, 0, 0, 0 )
! !- inverse of RGB attenuation lengths
zlogc = zc(ji,jj,0)
zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
zCze = zc(ji,jj,1)
zrdpsi = zc(ji,jj,2) ! 1/zdelpsi
zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm) ! gdepw/zze
!!st05 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze
! ! make sure zchl value is such that: 0.03 < zchl < 10.
zchl = MAX( 0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp ) )
! ! Convert chlorophyll value to attenuation coefficient
irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) ! look-up table index
! Green ! Blue
r1_LG = rkrgb(2,irgb) ; r1_LB = rkrgb(1,irgb)
!
! !- fluxes at jk+1 w-level
ze3t = e3t(ji,jj,jk,Kmm)
zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue
zzeT = ( zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - -
#if defined key_RK3
! !- RK3 : temperature trend at jk t-level
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
#else
! !- MLF : heat content trend due to Qsr flux (qsr_hc)
qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
#endif
zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue store at jk+1 w-level
zeT(ji,jj) = zzeT ! total - - -
END_2D
END DO
!
DO jk = nkG+1, nkB !* down to Blue extinction *! (< ~1300 m : B , IR+RG removed from calculation)
DO_2D( 0, 0, 0, 0 )
! !- inverse of RGB attenuation lengths
zlogc = zc(ji,jj,0)
zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
zCze = zc(ji,jj,1)
zrdpsi = zc(ji,jj,2) ! 1/zdelpsi
zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm) ! gdepw/zze
!!st05 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze
! ! make sure zchl value is such that: 0.03 < zchl < 10.
zchl = MAX( 0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp ) )
! ! Convert chlorophyll value to attenuation coefficient
irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) ! look-up table index
r1_LB = rkrgb(1,irgb) ! Blue
!
! !- fluxes at jk+1 w-level
ze3t = e3t(ji,jj,jk,Kmm)
zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Blue
zzeT = ( zzeB ) * wmask(ji,jj,jk+1) ! Total - -
#if defined key_RK3
! !- RK3 : temperature trend at jk t-level
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
#else
! !- MLF : heat content trend due to Qsr flux (qsr_hc)
qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
#endif
zeB(ji,jj) = zzeB ! Blue store at jk+1 w-level
zeT(ji,jj) = zzeT ! total - - -
END_2D
END DO
!
END SUBROUTINE qsr_RGBc
SUBROUTINE qsr_RGB( kt, Kmm, pts, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE qsr_RGB ***
!!
!! ** Purpose : Red-Green-Blue solar radiation with constant chlorophyll
!!
!! ** Method : The profile of the solar radiation within the ocean is defined
!! through 2 wavebands (rn_si0,rn_si1) or 1 (rn_si0,rn_abs) + 3 wavebands (RGB)
!! At the bottom, boudary condition for the radiation is no flux :
!! all heat which has not been absorbed in the above levels is put
!! in the last ocean level.
!! For each band, the computation is only done down to the level where
!! I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) .
!!
!! ** Action : - update ta with the penetrative solar radiation trend
!! - send trend for further diagnostics (l_trdtra=T)
!!
!! Reference : Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
!! Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
!!----------------------------------------------------------------------
INTEGER, INTENT(in ) :: kt, Kmm, Krhs ! ocean time-step and time-level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
!!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zze0, zzeR, zzeG, zzeB, zzeT ! - -
REAL(wp) :: zz0 , zz1 , ze3t ! - -
REAL(wp), DIMENSION(A2D(0)) :: ze0, zeR, zeG, zeB, zeT
!!----------------------------------------------------------------------
!
!
! !==============================================!
! !== R-G-B fluxes with constant chlorophyll ==!
! !======================********================!
!
! != surface light =!
!
zz0 = rn_abs ! Infrared absorption
zz1 = ( 1._wp - rn_abs ) / 3._wp ! surface equi-partition in R-G-B
!
DO_2D( 0, 0, 0, 0 ) ! surface light
ze0(ji,jj) = zz0 * qsr(ji,jj) ; zeR(ji,jj) = zz1 * qsr(ji,jj) ! IR ; Red
zeG(ji,jj) = zz1 * qsr(ji,jj) ; zeB(ji,jj) = zz1 * qsr(ji,jj) ! Green ; Blue
zeT(ji,jj) = qsr(ji,jj) ! Total
END_2D
!
! != interior light =!
!
DO jk = 1, nk0 !* near surface layers *! (< ~12 meters : IR + RGB )
DO_2D( 0, 0, 0, 0 )
ze3t = e3t(ji,jj,jk,Kmm)
zze0 = ze0(ji,jj) * EXP( - ze3t * r1_si0 ) ; zzeR = zeR(ji,jj) * EXP( - ze3t * r1_LR ) ! IR ; Red at jk+1 w-level
zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue - -
zzeT = ( zze0 + zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1) ! Total - -
!!st7-9 zzeT = ( zze0 + zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - -
#if defined key_RK3
! ! RK3 : temperature trend at jk t-level
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
#else
! ! MLF : heat content trend due to Qsr flux (qsr_hc)
qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
#endif
ze0(ji,jj) = zze0 ; zeR(ji,jj) = zzeR ! IR ; Red store at jk+1 w-level
zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue - - -
zeT(ji,jj) = zzeT ! total - - -
END_2D
!!stbug IF( jk == 1 ) THEN !* sea-ice *! store the 1st level attenuation coeff.
!!stbug WHERE( qsr(A2D(0)) /= 0._wp ) ; fraqsr_1lev(A2D(0)) = 1._wp - zeT(A2D(0)) / qsr(A2D(0))
!!stbug ELSEWHERE ; fraqsr_1lev(A2D(0)) = 1._wp
!!stbug END WHERE
!!stbug ENDIF
END DO
!
DO jk = nk0+1, nkR !* down to Red extinction *! (< ~71 meters : RGB , IR removed from calculation)
DO_2D( 0, 0, 0, 0 )
ze3t = e3t(ji,jj,jk,Kmm)
zzeR = zeR(ji,jj) * EXP( - ze3t * r1_LR ) ! Red at jk+1 w-level
zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue - -
zzeT = ( zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1) ! Total - -
!!st7-11 zzeT = ( zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - -
#if defined key_RK3
! ! RK3 : temperature trend at jk t-level
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
#else
! ! MLF : heat content trend due to Qsr flux (qsr_hc)
qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
#endif
zeR(ji,jj) = zzeR ! Red store at jk+1 w-level
zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue - - -
zeT(ji,jj) = zzeT ! total - - -
END_2D
END DO
!
DO jk = nkR+1, nkG !* down to Green extinction *! (< ~350 m : GB , IR+R removed from calculation)
DO_2D( 0, 0, 0, 0 )
ze3t = e3t(ji,jj,jk,Kmm)
zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue at jk+1 w-level
zzeT = ( zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - -
#if defined key_RK3
! ! RK3 : temperature trend at jk t-level
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
#else
! ! MLF : heat content trend due to Qsr flux (qsr_hc)
qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
#endif
zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue store at jk+1 w-level
zeT(ji,jj) = zzeT ! total - - -
END_2D
END DO
!
DO jk = nkG+1, nkB !* down to Blue extinction *! (< ~1300 m : B , IR+RG removed from calculation)
DO_2D( 0, 0, 0, 0 )
ze3t = e3t(ji,jj,jk,Kmm)
zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Blue at jk+1 w-level
zzeT = ( zzeB ) * wmask(ji,jj,jk+1) ! Total - -
#if defined key_RK3
! ! RK3 : temperature trend at jk t-level
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t
#else
! ! MLF : heat content trend due to Qsr flux (qsr_hc)
qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT )
#endif
zeB(ji,jj) = zzeB ! Blue store at jk+1 w-level
zeT(ji,jj) = zzeT ! total - - -
END_2D
END DO
!
END SUBROUTINE qsr_RGB
SUBROUTINE qsr_2BD( Kmm, pts, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE qsr_2BD ***
!!
!! ** Purpose : 2 bands (IR+visible) solar radiation with constant chlorophyll
!!
!! ** Method : The profile of the solar radiation within the ocean is defined
!! through 2 wavebands (rn_si0,rn_si1) a ratio rn_abs for IR absorbtion.
!! Considering the 2 wavebands case:
!! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
!! The temperature trend associated with the solar radiation penetration
!! is given by : zta = 1/e3t dk[ I ] / (rho0*Cp)
!! At the bottom, boudary condition for the radiation is no flux :
!! all heat which has not been absorbed in the above levels is put
!! in the last ocean level.
!! The computation is only done down to the level where
!! I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) .
!!
!! ** Action : - update ta with the penetrative solar radiation trend
!! - send trend for further diagnostics (l_trdtra=T)
!!
!! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
!! Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
!! Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
!!----------------------------------------------------------------------
INTEGER, INTENT(in ) :: Kmm, Krhs ! ocean time-step and time-level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
!!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zzatt ! - -
REAL(wp) :: zz0 , zz1 , ze3t ! - -
REAL(wp), DIMENSION(A2D(0)) :: zatt
!!----------------------------------------------------------------------
!
! !======================!
! !== 2-bands fluxes ==!
! !======================!
!
zz0 = rn_abs * r1_rho0_rcp ! surface equi-partition in 2-bands
zz1 = ( 1._wp - rn_abs ) * r1_rho0_rcp
!
zatt(A2D(0)) = r1_rho0_rcp !* surface value *!
!
DO_2D( 0, 0, 0, 0 )
zatt(ji,jj) = ( zz0 * EXP( -gdepw(ji,jj,1,Kmm)*r1_si0 ) + zz1 * EXP( -gdepw(ji,jj,1,Kmm)*r1_si1 ) )
END_2D
!
!!st IF(lwp) WRITE(numout,*) 'level = ', 1, ' qsr max = ' , MAXVAL(zatt)*rho0_rcp, ' W/m2', ' qsr min = ' , MINVAL(zatt)*rho0_rcp, ' W/m2'
!
DO jk = 1, nk0 !* near surface layers *! (< ~14 meters : IR + visible light )
DO_2D( 0, 0, 0, 0 )
ze3t = e3t(ji,jj,jk,Kmm) ! light attenuation at jk+1 w-level (divided by rho0_rcp)
zzatt = ( zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si0 ) &
& + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si1 ) ) * wmask(ji,jj,jk+1)
#if defined key_RK3
! ! RK3 : temperature trend at jk t-level
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) / ze3t
#else
! ! MLF : heat content trend due to Qsr flux (qsr_hc)
qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zatt(ji,jj) - zzatt )
#endif
zatt(ji,jj) = zzatt ! save for the next level computation
END_2D
!!stbug ! !* sea-ice *! store the 1st level attenuation coeff.
!!stbug IF( jk == 1 ) fraqsr_1lev(A2D(0)) = 1._wp - zatt(A2D(0)) * rho0_rcp
END DO
!!st IF(lwp) WRITE(numout,*) 'nk0+1= ', nk0+1, ' qsr max = ' , MAXVAL(zatt*qsr)*rho0_rcp, ' W/m2' , MAXVAL(zatt*qsr/e3t(:,:,nk0+1,Kmm)), ' K/s'
!
DO jk = nk0+1, nkV !* deeper layers *! (visible light only)
DO_2D( 0, 0, 0, 0 )
ze3t = e3t(ji,jj,jk,Kmm) ! light attenuation at jk+1 w-level (divided by rho0_rcp)
zzatt = ( zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si1 ) ) * wmask(ji,jj,jk+1)
#if defined key_RK3
! ! RK3 : temperature trend at jk t-level
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) / ze3t
#else
! ! MLF : heat content trend due to Qsr flux (qsr_hc)
qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zatt(ji,jj) - zzatt )
#endif
zatt(ji,jj) = zzatt ! save for the next level computation
END_2D
END DO
!
!!st IF(lwp) WRITE(numout,*) 'nkV+1= ', nkV+1, ' qsr max = ' , MAXVAL(zatt*qsr)*rho0_rcp, ' W/m2' , MAXVAL(zatt*qsr/e3t(:,:,nkV+1,Kmm)), ' K/s'
END SUBROUTINE qsr_2bd
FUNCTION qsr_ext_lev( pL, pfr ) RESULT( klev )
!!----------------------------------------------------------------------
!! *** ROUTINE trc_oce_ext_lev ***
!!
!! ** Purpose : compute the maximum level of light penetration
!!
!! ** Method : the function provides the level at which irradiance, I,
!! has a negligible effect on temperature.
!! T(n+1)-T(n) = ∆t dk[I] / ( rho0 Cp e3t_k )
!! I(k) has a negligible effect on temperature at level k if:
!! ∆t I(k) / ( rho0*Cp*e3t_k ) <= 1.e-15 °C
!! with I(z) = Qsr*pfr*EXP(-z/L), therefore :
!! z >= L * LOG( 1.e-15 * rho0*Cp*e3t_k / ( ∆t*Qsr*pfr ) )
!! with Qsr being the maximum normal surface irradiance at sea
!! level (~1000 W/m2).
!! # pL is the longest depth of extinction:
!! - pL = 23.00 m (2 bands case)
!! - pL = 48.24 m (3 bands case: blue waveband & 0.03 mg/m2 for the chlorophyll)
!! # pfr is the fraction of solar radiation which penetrates,
!! considering Qsr=1000 W/m2 and rn_abs = 0.58:
!! - Qsr*pfr0 = Qsr * rn_abs = 580 W/m2 (top absorbtion)
!! - Qsr*pfr1 = Qsr * (1-rn_abs) = 420 W/m2 (2 bands case)
!! - Qsr*pfr1 = Qsr * (1-rn_abs)/3 = 140 W/m2 (3 bands case & equi-partition)
!!
!!----------------------------------------------------------------------
INTEGER :: klev ! result: maximum level of light penetration
REAL(wp), INTENT(in) :: pL ! depth of extinction
REAL(wp), INTENT(in) :: pfr ! frac. solar radiation which penetrates
!
INTEGER :: jk ! dummy loop index
REAL(wp) :: zcoef ! local scalar
REAL(wp) :: zhext ! deepest depth until which light penetrates
REAL(wp) :: ze3t , zdw ! max( e3t_k ) and min( w-depth_k+1 )
REAL(wp) :: zprec = 10.e-15_wp ! required precision
REAL(wp) :: zQmax= 1000._wp ! maximum normal surface irradiance at sea level (W/m2)
!!----------------------------------------------------------------------
!
zQmax = 1000._wp ! maximum normal surface irradiance at sea level (W/m2)
!
zcoef = zprec * rho0_rcp / ( rDt * zQmax * pfr)
!
IF( ln_zco .OR. ln_zps ) THEN ! z- or zps coordinate (use 1D ref vertcial coordinate)
klev = jpkm1 ! Level of light extinction zco / zps
DO jk = jpkm1, 1, -1
zdw = gdepw_1d(jk+1) ! max w-depth at jk+1 level
ze3t = e3t_1d(jk ) ! minimum e3t at jk level
zhext = - pL * LOG( zcoef * ze3t ) ! extinction depth
IF( zdw >= zhext ) klev = jk ! last T-level reached by Qsr
END DO
ELSE ! s- or s-z- coordinate (use 3D vertical coordinate)
klev = jpkm1 ! Level of light extinction
DO jk = jpkm1, 1, -1 !
IF( SUM( tmask(:,:,jk) ) > 0 ) THEN ! ocean point at that level
zdw = MAXVAL( gdepw_0(:,:,jk+1) * wmask(:,:,jk) ) ! max w-depth at jk+1 level
ze3t = MINVAL( e3t_0(:,:,jk ) , mask=(wmask(:,:,jk+1)==1) ) ! minimum e3t at jk level
zhext = - pL * LOG( zcoef * ze3t ) ! extinction depth
IF( zdw >= zhext ) klev = jk ! last T-level reached by Qsr
ELSE ! only land point at level jk
klev = jk ! local domain sea-bed level
ENDIF
END DO
CALL mpp_max('tra_qsr', klev) ! needed for reproducibility !!st may be modified to avoid this comm.
! !!st use ssmask to remove the comm ?
ENDIF
!
!!st IF(lwp) WRITE(numout,*) ' level of e3t light extinction = ', klev, ' ref depth = ', gdepw_1d(klev+1), ' m'
END FUNCTION qsr_ext_lev
SUBROUTINE tra_qsr_init
!!----------------------------------------------------------------------
!! *** ROUTINE tra_qsr_init ***
......@@ -340,9 +776,9 @@ CONTAINS
!! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ios, irgb, ierror, ioptio ! local integer
REAL(wp) :: zz0, zc0 , zc1, zcoef ! local scalars
REAL(wp) :: zz1, zc2 , zc3, zchl ! - -
INTEGER :: ios, ierror, ioptio ! local integer
REAL(wp) :: zLB, zLG, zLR ! local scalar
REAL(wp) :: zVlp, zchl ! - -
!
CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files
TYPE(FLD_N) :: sn_chl ! informations about the chlorofyl field to be read
......@@ -353,12 +789,11 @@ CONTAINS
!
READ ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' )
!
READ ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
READ ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902)
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' )
IF(lwm) WRITE ( numond, namtra_qsr )
!
IF(lwp) THEN ! control print
IF(lwp) THEN !** control print **!
WRITE(numout,*)
WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
WRITE(numout,*) '~~~~~~~~~~~~'
......@@ -368,12 +803,12 @@ CONTAINS
WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio
WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta
WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs
WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0
WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1
WRITE(numout,*) ' RGB & 2 bands: shortess attenuation depth rn_si0 = ', rn_si0
WRITE(numout,*) ' 2 bands: longest attenuation depth rn_si1 = ', rn_si1
WRITE(numout,*)
ENDIF
!
ioptio = 0 ! Parameter control
ioptio = 0 !** Parameter control **!
IF( ln_qsr_rgb ) ioptio = ioptio + 1
IF( ln_qsr_2bd ) ioptio = ioptio + 1
IF( ln_qsr_bio ) ioptio = ioptio + 1
......@@ -386,21 +821,50 @@ CONTAINS
IF( ln_qsr_2bd ) nqsr = np_2BD
IF( ln_qsr_bio ) nqsr = np_BIO
!
! ! Initialisation
xsi0r = 1._wp / rn_si0
xsi1r = 1._wp / rn_si1
! !** Initialisation **!
!
! !== Infrared attenuation ==! (all schemes)
! !============================!
!
r1_si0 = 1._wp / rn_si0 ! inverse of infrared attenuation length
!
nk0 = qsr_ext_lev( rn_si0, rn_abs ) ! level of light extinction
!
IF(lwp) WRITE(numout,*) ' ==>>> Infrared light attenuation'
IF(lwp) WRITE(numout,*) ' level of infrared extinction = ', nk0, ' ref depth = ', gdepw_1d(nk0+1), ' m'
IF(lwp) WRITE(numout,*)
!
SELECT CASE( nqsr )
!
CASE( np_RGB , np_RGBc ) !== Red-Green-Blue light penetration ==!
CASE( np_RGBc, np_RGB ) !== Red-Green-Blue light attenuation ==! (Chl data or constant)
! !========================================!
!
IF(lwp) WRITE(numout,*) ' ==>>> R-G-B light penetration '
IF( nqsr == np_RGB ) THEN ; zchl = 0.05 ! constant Chl value
ELSE ; zchl = 0.03 ! minimum Chl value
ENDIF
zchl = MAX( 0.03_wp , MIN( zchl , 10._wp) ) ! NB. make sure that chosen value verifies: 0.03 < zchl < 10
nc_rgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) ! Convert Chl value to attenuation coefficient look-up table index
!
CALL trc_oce_rgb( rkrgb ) ! tabulated attenuation coef.
!
nksr = trc_oce_ext_lev( r_si2, 33._wp ) ! level of light extinction
zVlp = ( 1._wp - rn_abs ) / 3._wp ! visible light equi-partition
!
IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
! 1 / length ! attenuation length ! attenuation level
r1_LR = rkrgb(3,nc_rgb) ; zLR = 1._wp / r1_LR ; nkR = qsr_ext_lev( zLR, zVlp ) ! Red
r1_LG = rkrgb(2,nc_rgb) ; zLG = 1._wp / r1_LG ; nkG = qsr_ext_lev( zLG, zVlp ) ! Green
r1_LB = rkrgb(1,nc_rgb) ; zLB = 1._wp / r1_LB ; nkB = qsr_ext_lev( zLB, zVlp ) ! Blue
!
nkV = nkB ! maximum level of light penetration
!
IF( nqsr == np_RGB ) THEN
IF(lwp) WRITE(numout,*) ' ==>>> RGB: light attenuation with a constant Chlorophyll = ', zchl
ELSE
IF(lwp) WRITE(numout,*) ' ==>>> RGB: light attenuation using Chlorophyll data with min(Chl) = ', zchl
ENDIF
IF(lwp) WRITE(numout,*) ' level of Red extinction = ', nkR, ' ref depth = ', gdepw_1d(nkR+1), ' m'
IF(lwp) WRITE(numout,*) ' level of Green extinction = ', nkG, ' ref depth = ', gdepw_1d(nkG+1), ' m'
IF(lwp) WRITE(numout,*) ' level of Blue extinction = ', nkB, ' ref depth = ', gdepw_1d(nkB+1), ' m'
IF(lwp) WRITE(numout,*)
!
IF( nqsr == np_RGBc ) THEN ! Chl data : set sf_chl structure
IF(lwp) WRITE(numout,*) ' ==>>> Chlorophyll read in a file'
......@@ -408,22 +872,25 @@ CONTAINS
IF( ierror > 0 ) THEN
CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' ) ; RETURN
ENDIF
ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) )
ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) )
IF( sn_chl%ln_tint ) ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
! ! fill sf_chl with sn_chl and control print
CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init', &
CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init', &
& 'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
ENDIF
IF( nqsr == np_RGB ) THEN ! constant Chl
IF(lwp) WRITE(numout,*) ' ==>>> Constant Chlorophyll concentration = 0.05'
ENDIF
!
CASE( np_2BD ) !== 2 bands light penetration ==!
CASE( np_2BD ) !== 2 bands light attenuation (IR+ visible light) ==!
!
IF(lwp) WRITE(numout,*) ' ==>>> 2 bands light penetration'
!
nksr = trc_oce_ext_lev( rn_si1, 100._wp ) ! level of light extinction
IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
IF( lk_top ) CALL trc_oce_rgb( rkrgb ) ! tabulated attenuation coef.
!
r1_si1 = 1._wp / rn_si1 ! inverse of visible light attenuation
zVlp = ( 1._wp - rn_abs ) ! visible light partition
nkV = qsr_ext_lev( rn_si1, zVlp ) ! level of visible light extinction
!
IF(lwp) WRITE(numout,*) ' ==>>> 2 bands attenuation (Infrared + Visible light) '
IF(lwp) WRITE(numout,*) ' level of visible light extinction = ', nkV, ' ref depth = ', gdepw_1d(nkV+1), ' m'
IF(lwp) WRITE(numout,*)
!
CASE( np_BIO ) !== BIO light penetration ==!
!
......@@ -432,15 +899,19 @@ CONTAINS
!
CALL trc_oce_rgb( rkrgb ) ! tabulated attenuation coef.
!
nksr = trc_oce_ext_lev( r_si2, 33._wp ) ! level of light extinction
nkV = trc_oce_ext_lev( r_si2, 33._wp ) ! maximum level of light extinction
!
IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
IF(lwp) WRITE(numout,*) ' level of light extinction = ', nkV, ' ref depth = ', gdepw_1d(nkV+1), ' m'
!
END SELECT
!
qsr_hc(:,:,:) = 0._wp ! now qsr heat content set to zero where it will not be computed
nksr = nkV ! name of max level of light extinction used in traatf(_qco).F90
!
#if ! defined key_RK3
qsr_hc(:,:,:) = 0._wp ! MLF : now qsr heat content set to zero where it will not be computed
#endif
!
! 1st ocean level attenuation coefficient (used in sbcssm)
! ! Sea-ice : 1st ocean level attenuation coefficient (used in sbcssm)
IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
CALL iom_get( numror, jpdom_auto, 'fraqsr_1lev' , fraqsr_1lev )
ELSE
......
src/OCE/TRA/trasbc.F90
View file @ d1c6ae27
......@@ -38,14 +38,15 @@ MODULE trasbc
IMPLICIT NONE
PRIVATE
PUBLIC tra_sbc ! routine called by step.F90
PUBLIC tra_sbc ! routine called by step.F90
PUBLIC tra_sbc_RK3 ! routine called by stprk3_stg.F90
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: trasbc.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: trasbc.F90 15379 2021-10-15 09:05:45Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -221,6 +222,168 @@ CONTAINS
IF( ln_timing ) CALL timing_stop('tra_sbc')
!
END SUBROUTINE tra_sbc
SUBROUTINE tra_sbc_RK3 ( kt, Kmm, pts, Krhs, kstg )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_sbc_RK3 ***
!!
!! ** Purpose : Compute the tracer surface boundary condition trend of
!! (flux through the interface, concentration/dilution effect)
!! and add it to the general trend of tracer equations.
!!
!! ** Method : The (air+ice)-sea flux has two components:
!! (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface);
!! (2) Fwe , tracer carried with the water that is exchanged with air+ice.
!! The input forcing fields (emp, rnf, sfx) contain Fext+Fwe,
!! they are simply added to the tracer trend (ts(Krhs)).
!! In linear free surface case (ln_linssh=T), the volume of the
!! ocean does not change with the water exchanges at the (air+ice)-sea
!! interface. Therefore another term has to be added, to mimic the
!! concentration/dilution effect associated with water exchanges.
!!
!! ** Action : - Update ts(Krhs) with the surface boundary condition trend
!! - send trends to trdtra module for further diagnostics(l_trdtra=T)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt, Kmm, Krhs ! ocean time-step and time-level indices
INTEGER , INTENT(in ) :: kstg ! RK3 stage index
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer Eq.
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp) :: z1_rho0_e3t, zdep, ztim ! local scalar
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('tra_sbc_RK3')
!
IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile
!
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'tra_sbc_RK3 : TRAcer Surface Boundary Condition'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
ENDIF
!
IF( l_trdtra ) THEN !* Save ta and sa trends
ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
ENDIF
!
ENDIF
!
!!gm This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
IF( .NOT.ln_traqsr .AND. kstg == 1) THEN ! no solar radiation penetration
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
qns(ji,jj) = qns(ji,jj) + qsr(ji,jj) ! total heat flux in qns
qsr(ji,jj) = 0._wp ! qsr set to zero
END_2D
ENDIF
!----------------------------------------
! EMP, SFX and QNS effects
!----------------------------------------
! !== update tracer trend ==!
SELECT CASE( kstg )
!
CASE( 1 , 2 ) != stage 1 and 2 =! only in non linear ssh
!
IF( .NOT.ln_linssh ) THEN !* only heat and salt fluxes associated with mass fluxes
DO_2D( 0, 0, 0, 0 )
z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm)
pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) - emp(ji,jj)*pts(ji,jj,1,jp_tem,Kmm) * z1_rho0_e3t
pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) - emp(ji,jj)*pts(ji,jj,1,jp_sal,Kmm) * z1_rho0_e3t
END_2D
ENDIF
!
CASE( 3 )
!
IF( ln_linssh ) THEN !* linear free surface
DO_2D( 0, 0, 0, 0 )
z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm)
pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + ( r1_rcp * qns(ji,jj) & ! non solar heat flux
& + emp(ji,jj)*pts(ji,jj,1,jp_tem,Kmm) ) * z1_rho0_e3t ! add concentration/dilution effect due to constant volume cell
pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) + ( sfx(ji,jj) & ! salt flux due to freezing/melting
& + emp(ji,jj)*pts(ji,jj,1,jp_sal,Kmm) ) * z1_rho0_e3t ! add concentration/dilution effect due to constant volume cell
END_2D
IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile
IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) )
IF( iom_use('emp_x_sss') ) CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) )
ENDIF
ELSE
DO_2D( 0, 0, 0, 0 )
z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm)
pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + r1_rcp * qns(ji,jj) * z1_rho0_e3t
pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) + sfx(ji,jj) * z1_rho0_e3t
END_2D
ENDIF
END SELECT
!
!
!----------------------------------------
! River Runoff effects
!----------------------------------------
!
IF( ln_rnf ) THEN ! input of heat and salt due to river runoff
DO_2D( 0, 0, 0, 0 )
IF( rnf(ji,jj) /= 0._wp ) THEN
zdep = 1._wp / h_rnf(ji,jj)
DO jk = 1, nk_rnf(ji,jj)
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + rnf_tsc(ji,jj,jp_tem) * zdep
IF( ln_rnf_sal ) pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) + rnf_tsc(ji,jj,jp_sal) * zdep
END DO
ENDIF
END_2D
ENDIF
!
IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile
IF( iom_use('rnf_x_sst') ) CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) ) ! runoff term on sst
IF( iom_use('rnf_x_sss') ) CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) ) ! runoff term on sss
ENDIF
#if defined key_asminc
!
!----------------------------------------
! Assmilation effects
!----------------------------------------
!
IF( ln_sshinc .AND. kstg == 3 ) THEN ! input of heat and salt due to assimilation
!
IF( ln_linssh ) THEN
DO_2D( 0, 0, 0, 0 )
ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm)
pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + pts(ji,jj,1,jp_tem,Kmm) * ztim
pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) + pts(ji,jj,1,jp_sal,Kmm) * ztim
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
ztim = ssh_iau(ji,jj) / ( ht(ji,jj) + 1. - ssmask(ji, jj) )
pts(ji,jj,:,jp_tem,Krhs) = pts(ji,jj,:,jp_tem,Krhs) + pts(ji,jj,:,jp_tem,Kmm) * ztim
pts(ji,jj,:,jp_sal,Krhs) = pts(ji,jj,:,jp_sal,Krhs) + pts(ji,jj,:,jp_sal,Kmm) * ztim
END_2D
ENDIF
!
ENDIF
!
#endif
!
IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics
IF( ntile == 0 .OR. ntile == nijtile ) THEN
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt )
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds )
DEALLOCATE( ztrdt , ztrds )
ENDIF
ENDIF
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' sbc - Ta: ', mask1=tmask, &
& tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
!
IF( ln_timing ) CALL timing_stop('tra_sbc_RK3')
!
END SUBROUTINE tra_sbc_RK3
!!======================================================================
END MODULE trasbc
src/OCE/ZDF/zdfphy.F90
View file @ d1c6ae27
......@@ -209,7 +209,7 @@ CONTAINS
ioptio = 0
IF( ln_zdfcst ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_CST ; ENDIF
IF( ln_zdfric ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_RIC ; CALL zdf_ric_init ; ENDIF
IF( ln_zdftke ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_TKE ; CALL zdf_tke_init( Kmm ) ; ENDIF
IF( ln_zdftke ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_TKE ; CALL zdf_tke_init ; ENDIF
IF( ln_zdfgls ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_GLS ; CALL zdf_gls_init ; ENDIF
IF( ln_zdfosm ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_OSM ; CALL zdf_osm_init( Kmm ) ; ENDIF
!
......@@ -350,7 +350,7 @@ CONTAINS
IF( ln_zdfiwm ) CALL zdf_iwm( kt, Kmm, avm, avt, avs ) ! internal wave (de Lavergne et al 2017)
! !* Lateral boundary conditions (sign unchanged)
IF(nn_hls==1) THEN
IF(nn_hls==1) THEN ! if nn_hls==2 lbc_lnk done in stp routines
IF( l_zdfsh2 ) THEN
CALL lbc_lnk( 'zdfphy', avm_k, 'W', 1.0_wp , avt_k, 'W', 1.0_wp, &
& avm , 'W', 1.0_wp , avt , 'W', 1.0_wp , avs , 'W', 1.0_wp )
......
src/OCE/ZDF/zdftke.F90
View file @ d1c6ae27
......@@ -698,7 +698,7 @@ CONTAINS
END SUBROUTINE tke_avn
SUBROUTINE zdf_tke_init( Kmm )
SUBROUTINE zdf_tke_init
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_tke_init ***
!!
......@@ -714,7 +714,6 @@ CONTAINS
!!----------------------------------------------------------------------
USE zdf_oce , ONLY : ln_zdfiwm ! Internal Wave Mixing flag
!!
INTEGER, INTENT(in) :: Kmm ! time level index
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ios
!!
......
src/OCE/module_example.F90
View file @ d1c6ae27
......@@ -54,7 +54,7 @@ MODULE exampl
# include "exampl_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: module_example.F90 14842 2021-05-11 13:17:26Z acc $
!! $Id: module_example.F90 14941 2021-06-03 11:42:27Z acc $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -144,7 +144,7 @@ CONTAINS
!! ** Action : ...
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk, jit ! dummy loop indices
INTEGER :: ios ! Local integer output status for namelist read
INTEGER :: ios ! Local integer output status for namelist read
!!
NAMELIST/namexa/ exa_v1, exa_v2, nexa_0, sn_ex
!!----------------------------------------------------------------------
......
src/OCE/nemogcm.F90
View file @ d1c6ae27
......@@ -65,7 +65,11 @@ MODULE nemogcm
USE ice_domain_size, only: nx_global, ny_global
#endif
#if defined key_qco || defined key_linssh
# if defined key_RK3
USE stprk3
# else
USE stpmlf ! NEMO time-stepping (stp_MLF routine)
# endif
#else
USE step ! NEMO time-stepping (stp routine)
#endif
......@@ -91,7 +95,7 @@ MODULE nemogcm
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: nemogcm.F90 15267 2021-09-17 09:04:34Z smasson $
!! $Id: nemogcm.F90 15532 2021-11-24 11:47:32Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -121,7 +125,11 @@ CONTAINS
CALL nemo_init !== Initialisations ==!
! !-----------------------!
#if defined key_agrif
Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices
# if defined key_RK3
Kbb_a = Nbb; Kmm_a = Nbb; Krhs_a = Nrhs ! RK3: agrif_oce module copies of time level indices
# else
Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! MLF: agrif_oce module copies of time level indices
# endif
CALL Agrif_Declare_Var ! " " " " " DYN/TRA
# if defined key_top
CALL Agrif_Declare_Var_top ! " " " " " TOP
......@@ -146,14 +154,22 @@ CONTAINS
CALL Agrif_Regrid()
!
! Recursive update from highest nested level to lowest:
Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices
# if defined key_RK3
Kbb_a = Nbb; Kmm_a = Nbb; Krhs_a = Nrhs ! RK3: agrif_oce module copies of time level indices
# else
Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! MLF: agrif_oce module copies of time level indices
# endif
CALL Agrif_step_child_adj(Agrif_Update_All)
CALL Agrif_step_child_adj(Agrif_Check_parent_bat)
!
DO WHILE( istp <= nitend .AND. nstop == 0 )
!
# if defined key_qco || defined key_linssh
# if defined key_RK3
CALL stp_RK3
# else
CALL stp_MLF
# endif
# else
CALL stp
# endif
......@@ -174,7 +190,11 @@ CONTAINS
ENDIF
!
# if defined key_qco || defined key_linssh
# if defined key_RK3
CALL stp_RK3( istp )
# else
CALL stp_MLF( istp )
# endif
# else
CALL stp ( istp )
# endif
......@@ -391,7 +411,11 @@ CONTAINS
! Initialise time level indices
Nbb = 1 ; Nnn = 2 ; Naa = 3 ; Nrhs = Naa
#if defined key_agrif
Kbb_a = Nbb ; Kmm_a = Nnn ; Krhs_a = Nrhs ! agrif_oce module copies of time level indices
# if defined key_RK3
Kbb_a = Nbb ; Kmm_a = Nbb ; Krhs_a = Nrhs ! RK3: agrif_oce module copies of time level indices
# else
Kbb_a = Nbb ; Kmm_a = Nnn ; Krhs_a = Nrhs ! MLF: agrif_oce module copies of time level indices
# endif
#endif
! !-------------------------------!
! ! NEMO general initialization !
......@@ -470,7 +494,11 @@ CONTAINS
CALL isf_init( Nbb, Nnn, Naa )
#if defined key_top
! ! Passive tracers
# if defined key_RK3
CALL trc_init( Nbb, Nbb, Naa )
# else
CALL trc_init( Nbb, Nnn, Naa )
# endif
#endif
IF( l_ldfslp ) CALL ldf_slp_init ! slope of lateral mixing
......
src/OCE/oce.F90
View file @ d1c6ae27
......@@ -20,8 +20,8 @@ MODULE oce
!! dynamics and tracer fields
!! --------------------------
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uu , vv !: horizontal velocities [m/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ww !: vertical velocity [m/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:), TARGET :: uu , vv !: horizontal velocities [m/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) , TARGET :: ww !: vertical velocity [m/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wi !: vertical vel. (adaptive-implicit) [m/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv !: horizontal divergence [s-1]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: ts !: 4D T-S fields [Celsius,psu]
......@@ -45,10 +45,13 @@ MODULE oce
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_e !: external v-depth
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hur_e !: inverse of u-depth
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hvr_e !: inverse of v-depth
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv, vn_adv !: Advective barotropic fluxes
#if ! defined key_RK3
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_b , vb2_b !: Half step fluxes (ln_bt_fw=T)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_bf , vn_bf !: Asselin filtered half step fluxes (ln_bt_fw=T)
#endif
#if defined key_agrif
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_i_b, vb2_i_b !: Half step time integrated fluxes
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_i_b, vb2_i_b !: agrif time integrated fluxes
#endif
!! interpolated gradient (only used in zps case)
......@@ -76,7 +79,7 @@ MODULE oce
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: oce.F90 15556 2021-11-29 15:23:06Z jchanut $
!! $Id: oce.F90 14381 2021-02-03 12:36:25Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -85,37 +88,48 @@ CONTAINS
!!----------------------------------------------------------------------
!! *** FUNCTION oce_alloc ***
!!----------------------------------------------------------------------
INTEGER :: ierr(6)
INTEGER :: ii
INTEGER, DIMENSION(7) :: ierr
!!----------------------------------------------------------------------
!
ierr(:) = 0
ii = 0 ; ierr(:) = 0
!
ii = ii+1
ALLOCATE( uu (jpi,jpj,jpk,jpt) , vv (jpi,jpj,jpk,jpt) , &
& ww (jpi,jpj,jpk) , hdiv(jpi,jpj,jpk) , &
& ts (jpi,jpj,jpk,jpts,jpt) , &
& rab_b(jpi,jpj,jpk,jpts) , rab_n(jpi,jpj,jpk,jpts) , &
& rn2b (jpi,jpj,jpk) , rn2 (jpi,jpj,jpk) , &
& rhd (jpi,jpj,jpk) , rhop (jpi,jpj,jpk) , STAT=ierr(1) )
& rhd (jpi,jpj,jpk) , rhop (jpi,jpj,jpk) , STAT=ierr(ii) )
!
ii = ii+1
ALLOCATE( ssh (jpi,jpj,jpt) , uu_b(jpi,jpj,jpt) , vv_b(jpi,jpj,jpt) , &
& ssh_frc(jpi,jpj) , &
& gtsu(jpi,jpj,jpts) , gtsv(jpi,jpj,jpts) , &
& gru (jpi,jpj) , grv (jpi,jpj) , &
& gtui(jpi,jpj,jpts) , gtvi(jpi,jpj,jpts) , &
& grui(jpi,jpj) , grvi(jpi,jpj) , &
& riceload(jpi,jpj) , STAT=ierr(2) )
& riceload(jpi,jpj) , STAT=ierr(ii) )
!
ALLOCATE( fraqsr_1lev(jpi,jpj) , STAT=ierr(3) )
ii = ii+1
ALLOCATE( fraqsr_1lev(jpi,jpj) , STAT=ierr(ii) )
!
ii = ii+1
ALLOCATE( ssha_e(jpi,jpj), sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), &
& ua_e(jpi,jpj), un_e(jpi,jpj), ub_e(jpi,jpj), ubb_e(jpi,jpj), &
& va_e(jpi,jpj), vn_e(jpi,jpj), vb_e(jpi,jpj), vbb_e(jpi,jpj), &
& hu_e(jpi,jpj), hur_e(jpi,jpj), hv_e(jpi,jpj), hvr_e(jpi,jpj), STAT=ierr(4) )
& hu_e(jpi,jpj), hur_e(jpi,jpj), hv_e(jpi,jpj), hvr_e(jpi,jpj), &
& un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(ii) )
!
ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_bf(jpi,jpj), vn_bf(jpi,jpj) , STAT=ierr(5) )
#if ! defined key_RK3
ii = ii+1 ! MLF: arrays related to Asselin filter + ???
ALLOCATE( un_bf(jpi,jpj), vn_bf(jpi,jpj), ub2_b(jpi,jpj), vb2_b(jpi,jpj) , STAT=ierr(ii) )
#endif
#if defined key_agrif
ALLOCATE( ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , STAT=ierr(6) )
ii = ii+1 ! AGRIF: ???
ALLOCATE( ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , STAT=ierr(ii) )
#endif
!
!
oce_alloc = MAXVAL( ierr )
IF( oce_alloc /= 0 ) CALL ctl_stop( 'STOP', 'oce_alloc: failed to allocate arrays' )
!
......
src/OCE/stp2d.F90 0 → 100644
View file @ d1c6ae27
MODULE stp2d
!!======================================================================
!! *** MODULE stp2d ***
!! Time-stepping : manager of the ocean, tracer and ice time stepping
!! using a 3rd order Rung Kuta with fixed or quasi-eulerian coordinate
!!======================================================================
!! History : 4.5 ! 2021-01 (S. Techene, G. Madec, N. Ducousso, F. Lemarie) Original code
!! NEMO
!!----------------------------------------------------------------------
#if defined key_qco || defined key_linssh
!!----------------------------------------------------------------------
!! 'key_qco' Quasi-Eulerian vertical coordinate
!! OR
!! 'key_linssh Fixed in time vertical coordinate
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! stp_2D : RK3 case
!!----------------------------------------------------------------------
USE step_oce ! time stepping used modules
USE domqco ! quasi-eulerian coordinate (dom_qco_r3c routine)
USE dynspg_ts ! 2D mode integration
USE sshwzv ! vertical speed
USE sbc_ice , ONLY : snwice_mass, snwice_mass_b
USE sbcapr ! surface boundary condition: atmospheric pressure
USE sbcwave, ONLY : bhd_wave
#if defined key_agrif
USE agrif_oce_interp
USE agrif_oce_sponge
#endif
PRIVATE
PUBLIC stp_2D ! called by nemogcm.F90
REAL (wp) :: r1_2 = 0.5_wp
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE stp_2D( kt, Kbb, Kmm, Kaa, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE stp_2D ***
!!
!! ** Purpose : - Compute sea-surface height and barotropic velocity at Kaa
!! in single 1st RK3.
!!
!! ** Method : -1- Compute the 3D to 2D forcing
!! * Momentum (Ue,Ve)_rhs :
!! 3D to 2D dynamics, i.e. the vertical sum of :
!! - Hor. adv. : KEG + RVO in vector form
!! : ADV_h + MET in flux form
!! - LDF Lateral mixing
!! - HPG Hor. pressure gradient
!! External forcings
!! - baroclinic drag
!! - wind
!! - atmospheric pressure
!! - snow+ice load
!! - surface wave load
!! * ssh (sshe_rhs) :
!! Net column average freshwater flux
!!
!! -2- Solve the external mode Eqs. using sub-time step
!! by a call to dyn_spg_ts (will be renamed dyn_2D or stp_2D)
!!
!! ** action : ssh : N+1 sea surface height (Kaa=N+1)
!! (uu_b,vv_b) : N+1 barotropic velocity
!! (un_adv,vn_adv): barotropic transport from N to N+1
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt, Kbb, Kmm, Kaa, Krhs ! ocean time-step and time-level indices
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zg_2, zintp, zgrho0r, zld, zztmp ! local scalars
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpice ! 2D workspace
!! ---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('stp_2D')
!
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'stp_2D : barotropic field in single first '
IF(lwp) WRITE(numout,*) '~~~~~~'
ENDIF
!
IF( ln_linssh ) THEN !== Compute ww(:,:,1) ==! (needed for momentum advection)
!!gm only in Flux Form, Vector Form dzU_z=0 assumed to be zero
!!gm ww(k=1) = div_h(uu_b) ==> modif dans dynadv <<<=== TO BE DONE
ENDIF
ALLOCATE( sshe_rhs(jpi,jpj) , Ue_rhs(jpi,jpj) , Ve_rhs(jpi,jpj) , CdU_u(jpi,jpj) , CdU_v(jpi,jpj) )
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RHS of barotropic momentum Eq.
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! !======================================!
! !== Dynamics 2D RHS from 3D trends ==! (HADV + LDF + HPG) (No Coriolis trend)
! !======================================!
uu(:,:,:,Krhs) = 0._wp ! set dynamics trends to zero
vv(:,:,:,Krhs) = 0._wp
!
! !* compute advection + coriolis *!
!
CALL ssh_nxt( kt, Kbb, Kbb, ssh, Kaa )
!
IF( .NOT.lk_linssh ) THEN
DO_2D_OVR( 1, nn_hls, 1, nn_hls ) ! loop bounds limited by ssh definition in ssh_nxt
r3t(ji,jj,Kaa) = ssh(ji,jj,Kaa) * r1_ht_0(ji,jj) ! "after" ssh/h_0 ratio guess at t-column at Kaa (n+1)
END_2D
ENDIF
!
CALL wzv ( kt, Kbb, Kbb, Kaa , uu(:,:,:,Kbb), vv(:,:,:,Kbb), ww ) ! ww guess at Kbb (n)
!
CALL dyn_adv( kt, Kbb, Kbb , uu, vv, Krhs) !- vector form KEG+ZAD
! !- flux form ADV
CALL dyn_vor( kt, Kbb, uu, vv, Krhs ) !- vector form COR+RVO
! !- flux form COR+MET
!
! !* lateral viscosity *!
CALL dyn_ldf( kt, Kbb, Kbb, uu, vv, Krhs )
#if defined key_agrif
IF(.NOT. Agrif_Root() ) THEN !* AGRIF: sponge *!
CALL Agrif_Sponge_dyn
ENDIF
#endif
!
! !* hydrostatic pressure gradient *!
CALL eos ( ts , Kbb, rhd ) ! in situ density anomaly at Kbb
CALL dyn_hpg( kt , Kbb , uu, vv, Krhs ) ! horizontal gradient of Hydrostatic pressure
!
! !* vertical averaging *!
Ue_rhs(:,:) = SUM( e3u_0(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu_0(:,:)
Ve_rhs(:,:) = SUM( e3v_0(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv_0(:,:)
! !===========================!
! !== external 2D forcing ==!
! !===========================!
!
! !* baroclinic drag forcing *! (also provide the barotropic drag coeff.)
!
CALL dyn_drg_init( Kbb, Kbb, uu, vv, uu_b, vv_b, Ue_rhs, Ve_rhs, CdU_u, CdU_v )
!
! !* wind forcing *!
IF( ln_bt_fw ) THEN
DO_2D( 0, 0, 0, 0 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kbb)
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kbb)
END_2D
ELSE
zztmp = r1_rho0 * r1_2
DO_2D( 0, 0, 0, 0 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kbb)
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kbb)
END_2D
ENDIF
!
! !* atmospheric pressure forcing *!
IF( ln_apr_dyn ) THEN
IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
DO_2D( 0, 0, 0, 0 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
END_2D
ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
zztmp = grav * r1_2
DO_2D( 0, 0, 0, 0 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
END_2D
ENDIF
ENDIF
!
! !* snow+ice load *! (embedded sea ice)
IF( ln_ice_embd ) THEN
ALLOCATE( zpice(jpi,jpj) )
zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
zgrho0r = - grav * r1_rho0
zpice(:,:) = ( zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:) ) * zgrho0r
DO_2D( 0, 0, 0, 0 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
END_2D
DEALLOCATE( zpice )
ENDIF
!
! !* surface wave load *! (Bernoulli head)
!
IF( ln_wave .AND. ln_bern_srfc ) THEN
DO_2D( 0, 0, 0, 0 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj)
END_2D
ENDIF
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RHS of see surface height Eq.
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
! !== Net water flux forcing ==! (applied to a water column)
!
IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
sshe_rhs(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) )
ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
zztmp = r1_rho0 * r1_2
sshe_rhs(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) &
& - rnf(:,:) - rnf_b(:,:) &
& - fwfisf_cav(:,:) - fwfisf_cav_b(:,:) &
& - fwfisf_par(:,:) - fwfisf_par_b(:,:) )
ENDIF
!
! !== Stokes drift divergence ==! (if exist)
!
IF( ln_sdw ) sshe_rhs(:,:) = sshe_rhs(:,:) + div_sd(:,:)
!
!
! !== ice sheet coupling ==!
!
IF( ln_isf .AND. ln_isfcpl ) THEN
IF( ln_rstart .AND. kt == nit000 ) sshe_rhs(:,:) = sshe_rhs(:,:) + risfcpl_ssh(:,:)
IF( ln_isfcpl_cons ) sshe_rhs(:,:) = sshe_rhs(:,:) + risfcpl_cons_ssh(:,:)
ENDIF
!
#if defined key_asminc
! !== Add the IAU weighted SSH increment ==!
!
IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) sshe_rhs(:,:) = sshe_rhs(:,:) - ssh_iau(:,:)
#endif
!
#if defined key_agrif
! !== AGRIF : fill boundary data arrays (on both )
IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
#endif
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Compute ssh and (uu_b,vv_b) at N+1 (Kaa)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! using a split-explicit time integration in forward mode
! ( ABM3-AM4 time-integration Shchepetkin et al. OM2005) with temporal diffusion (Demange et al. JCP2019) )
CALL dyn_spg_ts( kt, Kbb, Kbb, Krhs, uu, vv, ssh, uu_b, vv_b, Kaa ) ! time-splitting
DEALLOCATE( sshe_rhs , Ue_rhs , Ve_rhs , CdU_u , CdU_v )
!!gm this is useless I guess : RK3, done in each stages
!
! IF( ln_dynspg_ts ) THEN ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Krhs)
! ! as well as vertical scale factors and vertical velocity need to be updated
! CALL div_hor ( kstp, Kbb, Kmm ) ! Horizontal divergence (2nd call in time-split case)
! IF(.NOT.lk_linssh) CALL dom_qco_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa), r3f(:,:) ) ! update ssh/h_0 ratio at t,u,v,f pts
! ENDIF
!!gm
!
IF( ln_timing ) CALL timing_stop('stp_2D')
!
END SUBROUTINE stp_2D
#else
!!----------------------------------------------------------------------
!! default option EMPTY MODULE qco not activated
!!----------------------------------------------------------------------
#endif
!!======================================================================
END MODULE stp2d
src/OCE/stpmlf.F90
View file @ d1c6ae27
......@@ -34,7 +34,8 @@ MODULE stpmlf
!! 4.1 ! 2019-08 (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
!! 4.x ! 2020-08 (S. Techene, G. Madec) quasi eulerian coordinate time stepping
!!----------------------------------------------------------------------
#if defined key_qco || defined key_linssh
#if ! defined key_RK3
# if defined key_qco || defined key_linssh
!!----------------------------------------------------------------------
!! 'key_qco' Quasi-Eulerian vertical coordinate
!! OR
......@@ -91,8 +92,8 @@ CONTAINS
!! -7- Compute the diagnostics variables (rd,N2, hdiv,w)
!! -8- Outputs and diagnostics
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk, jtile ! dummy loop indice
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zgdept
INTEGER :: ji, jj, jk, jn, jtile ! dummy loop indice
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zgdept
!! ---------------------------------------------------------------------
#if defined key_agrif
IF( nstop > 0 ) RETURN ! avoid to go further if an error was detected during previous time step (child grid)
......@@ -281,7 +282,6 @@ CONTAINS
IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning
ENDIF
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! cool skin
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
......@@ -316,11 +316,17 @@ CONTAINS
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Active tracers
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
ts(:,:,:,:,Nrhs) = 0._wp ! set tracer trends to zero
IF( ln_tile ) CALL dom_tile_start ! [tiling] TRA tiling loop (1)
IF( ln_tile ) CALL dom_tile_start ! [tiling] TRA tiling loop (1)
DO jtile = 1, nijtile
IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
DO jn = 1, jpts
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
ts(ji,jj,jk,jn,Nrhs) = 0._wp ! set tracer trends to zero
END_3D
END DO
IF( lk_asminc .AND. ln_asmiau .AND. &
& ln_trainc ) CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs ) ! apply tracer assimilation increment
......@@ -473,7 +479,6 @@ CONTAINS
!! ** Action : puu(Kmm),pvv(Kmm) updated now horizontal velocity (ln_bt_fw=F)
!! puu(Kaa),pvv(Kaa) after horizontal velocity
!!----------------------------------------------------------------------
USE dynspg_ts, ONLY : un_adv, vn_adv ! updated Kmm barotropic transport
!!
INTEGER , INTENT(in ) :: Kmm, Kaa ! before and after time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities
......@@ -542,7 +547,7 @@ CONTAINS
! Update after tracer and velocity on domain lateral boundaries
!
# if defined key_agrif
CALL Agrif_tra !* AGRIF zoom boundaries
CALL Agrif_tra( kt ) !* AGRIF zoom boundaries
CALL Agrif_dyn( kt )
# endif
! ! local domain boundaries (T-point, unchanged sign)
......@@ -566,10 +571,11 @@ CONTAINS
!
END SUBROUTINE finalize_lbc
#else
# else
!!----------------------------------------------------------------------
!! default option EMPTY MODULE qco not activated
!!----------------------------------------------------------------------
# endif
#endif
!!======================================================================
......
src/OCE/stprk3.F90 0 → 100644
View file @ d1c6ae27
MODULE stprk3
!!======================================================================
!! *** MODULE stpRK3 ***
!! Time-stepping : manager of the ocean, tracer and ice time stepping
!! using a 3rd order Rung Kuta with fixed or quasi-eulerian coordinate
!!======================================================================
!! History : 4.5 ! 2021-01 (S. Techene, G. Madec, N. Ducousso, F. Lemarie) Original code
!! NEMO
!!----------------------------------------------------------------------
#if defined key_qco || defined key_linssh
!!----------------------------------------------------------------------
!! 'key_qco' Quasi-Eulerian vertical coordinate
!! OR
!! 'key_linssh Fixed in time vertical coordinate
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! stp_RK3 : NEMO 3rd order Runge-Kutta time-stepping
!!----------------------------------------------------------------------
USE step_oce ! time stepping used modules
USE domqco ! quasi-eulerian coordinate (dom_qco_r3c routine)
USE stprk3_stg ! RK3 stages
USE stp2d ! external mode solver
USE trd_oce ! trends: ocean variables
USE diaptr
USE ldftra
IMPLICIT NONE
PRIVATE
PUBLIC stp_RK3 ! called by nemogcm.F90
! !** time level indices **!
INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs !: used by nemo_init
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
#if defined key_agrif
RECURSIVE SUBROUTINE stp_RK3( )
INTEGER :: kstp ! ocean time-step index
#else
SUBROUTINE stp_RK3( kstp )
INTEGER, INTENT(in) :: kstp ! ocean time-step index
#endif
!!----------------------------------------------------------------------
!! *** ROUTINE stp_RK3 ***
!!
!! ** Purpose : - Time stepping of OCE (momentum and active tracer Eqs.) (RK3)
!! - Time stepping of SI3 (dynamic and thermodynamic Eqs.) (FBS)
!! - Time stepping of TRC (passive tracer Eqs.)
!!
!! ** Method : -1- Update forcings and data
!! -2- Update ocean physics
!! -3- Compute the after (Naa) ssh and velocity
!! -4- diagnostics and output at Now (Nnn)
!! -4- Compute the after (Naa) T-S
!! -5- Update now
!! -6- Update the horizontal velocity
!! -7- Compute the diagnostics variables (rd,N2, hdiv,w)
!! -8- Outputs and diagnostics
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk, jtile ! dummy loop indice
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zgdept
!! ---------------------------------------------------------------------
#if defined key_agrif
IF( nstop > 0 ) RETURN ! avoid to go further if an error was detected during previous time step (child grid)
kstp = nit000 + Agrif_Nb_Step()
Kbb_a = Nbb ; Kmm_a = Nbb ; Krhs_a = Nrhs ! agrif_oce module copies of time level indices
IF( lk_agrif_debug ) THEN
IF( Agrif_Root() .AND. lwp) WRITE(*,*) '---'
IF(lwp) WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
ENDIF
IF( kstp == nit000 + 1 ) lk_agrif_fstep = .FALSE.
# if defined key_xios
IF( Agrif_Nbstepint() == 0 ) CALL iom_swap( cxios_context )
# endif
#endif
!
IF( ln_timing ) CALL timing_start('stp_RK3')
!
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! update I/O and calendar
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
IF( kstp == nit000 ) THEN ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
CALL iom_init( cxios_context, ld_closedef=.FALSE. ) ! for model grid (including possible AGRIF zoom)
IF( lk_diamlr ) CALL dia_mlr_iom_init ! with additional setup for multiple-linear-regression analysis
CALL dia_ptr_init ! called here since it uses iom_use
CALL iom_init_closedef
IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid
ENDIF
IF( kstp == nitrst .AND. lwxios ) THEN
CALL iom_swap( cw_ocerst_cxt )
CALL iom_init_closedef( cw_ocerst_cxt )
CALL iom_setkt( kstp - nit000 + 1, cw_ocerst_cxt )
#if defined key_top
CALL iom_swap( cw_toprst_cxt )
CALL iom_init_closedef( cw_toprst_cxt )
CALL iom_setkt( kstp - nit000 + 1, cw_toprst_cxt )
#endif
ENDIF
#if defined key_si3
IF( kstp + nn_fsbc - 1 == nitrst .AND. lwxios ) THEN
CALL iom_swap( cw_icerst_cxt )
CALL iom_init_closedef( cw_icerst_cxt )
CALL iom_setkt( kstp - nit000 + 1, cw_icerst_cxt )
ENDIF
#endif
IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init)
CALL iom_setkt( kstp - nit000 + 1, cxios_context ) ! tell IOM we are at time step kstp
IF( ln_crs ) CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" ) ! tell IOM we are at time step kstp
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( ln_tide ) CALL tide_update( kstp ) ! update tide potential
IF( ln_apr_dyn ) CALL sbc_apr ( kstp ) ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
IF( ln_bdy ) CALL bdy_dta ( kstp, Nbb ) ! update dynamic & tracer data at open boundaries
IF( ln_isf ) CALL isf_stp ( kstp, Nbb ) ! update iceshelf geometry
CALL sbc ( kstp, Nbb, Nbb ) ! Sea Boundary Condition (including sea-ice)
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Update stochastic parameters and random T/S fluctuations
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
IF( ln_sto_eos ) CALL sto_par( kstp ) ! Stochastic parameters
IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn) ) ! Random T/S fluctuations
!!gm ocean physic computed at stage 3 ?
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Ocean physics update
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! THERMODYNAMICS
!!gm only Before is needed for bn2 and rab (except at stage 3 for n2)
!!gm issue with Nnn used in rad(Nbb)
CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nbb ) ! before local thermal/haline expension ratio at T-points
!!st CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn ) ! now local thermal/haline expension ratio at T-points
CALL bn2 ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nbb ) ! before Brunt-Vaisala frequency
!!st CALL bn2 ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency
!!gm
rab_n = rab_b
rn2 = rn2b
!!gm sh2 computed at the end of the time-step
!!gm or call zdf_phy at the end !
! VERTICAL PHYSICS
!!st CALL zdf_phy( kstp, Nbb, Nnn, Nrhs ) ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
CALL zdf_phy( kstp, Nbb, Nbb, Nrhs ) ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
!!gm gdep
! LATERAL PHYSICS
!
IF( ln_zps .OR. l_ldfslp ) CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) ) ! before in situ density
IF( ln_zps .AND. .NOT. ln_isfcav) &
& CALL zps_hde ( kstp, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient
& rhd, gru , grv ) ! of t, s, rd at the last ocean level
IF( ln_zps .AND. ln_isfcav) &
& CALL zps_hde_isf( kstp, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level
IF( l_ldfslp ) THEN ! slope of lateral mixing
IF( ln_traldf_triad ) THEN
CALL ldf_slp_triad( kstp, Nbb, Nbb ) ! before slope for triad operator
ELSE
CALL ldf_slp ( kstp, rhd, rn2b, Nbb, Nbb ) ! before slope for standard operator
ENDIF
ENDIF
! ! eddy diffusivity coeff.
IF( l_ldftra_time .OR. l_ldfeiv_time ) CALL ldf_tra( kstp, Nbb, Nbb ) ! and/or eiv coeff.
IF( l_ldfdyn_time ) CALL ldf_dyn( kstp, Nbb ) ! eddy viscosity coeff.
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RK3 : single first external mode computation
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
CALL stp_2D( kstp, Nbb, Nbb, Naa, Nrhs ) ! out: ssh, (uu_b,vv_b) at Naa and (un_adv,vn_adv) between Nbb and Naa
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RK3 time integration
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
CALL rk3_dia( kstp, 0 ) ! Diagnostics switched off for stage 1 & 2
!
! Stage 1 :
CALL stp_RK3_stg( 1, kstp, Nbb, Nbb, Nrhs, Naa )
!
Nrhs = Nnn ; Nnn = Naa ; Naa = Nrhs ! Swap: Nbb unchanged, Nnn <==> Naa
!
! Stage 2 :
CALL stp_RK3_stg( 2, kstp, Nbb, Nnn, Nrhs, Naa )
!
Nrhs = Nnn ; Nnn = Naa ; Naa = Nrhs ! Swap: Nbb unchanged, Nnn <==> Naa
!
! Stage 3 :
CALL rk3_dia( kstp, 1 ) ! Diagnostics switched on for stage 3
!
CALL stp_RK3_stg( 3, kstp, Nbb, Nnn, Nrhs, Naa )
!
Nrhs = Nbb ; Nbb = Naa ; Naa = Nrhs ! Swap: Nnn unchanged, Nbb <==> Naa
! linear extrapolation of ssh to compute ww at the beginning of the next time-step
! ssh(n+1) = 2*ssh(n) - ssh(n-1)
ssh(:,:,Naa) = 2*ssh(:,:,Nbb) - ssh(:,:,Naa)
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! diagnostics and outputs
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!==>>> at Nbb no more Nnn
IF( ln_diacfl ) CALL dia_cfl ( kstp, Nbb ) ! Courant number diagnostics
CALL dia_hth ( kstp, Nbb ) ! Thermocline depth (20 degres isotherm depth)
IF( ln_diadct ) CALL dia_dct ( kstp, Nbb ) ! Transports
!!st CALL dia_ar5 ( kstp, Nbb ) ! ar5 diag
CALL dia_ptr ( kstp, Nbb ) ! Poleward adv/ldf TRansports diagnostics
CALL dia_wri ( kstp, Nbb ) ! ocean model: outputs
IF( ln_crs ) CALL crs_fld ( kstp, Nbb ) ! ocean model: online field coarsening & output
IF( lk_diadetide ) CALL dia_detide( kstp ) ! Weights computation for daily detiding of model diagnostics
IF( lk_diamlr ) CALL dia_mlr ! Update time used in multiple-linear-regression analysis
!!====>>>> to be modified for RK3
! IF( ln_floats ) CALL flo_stp ( kstp, Nbb, Nnn ) ! drifting Floats
! IF( ln_diahsb ) CALL dia_hsb ( kstp, Nbb, Nnn ) ! - ML - global conservation diagnostics
!!st
IF( ln_diahsb ) CALL dia_hsb ( kstp, Nbb, Nbb ) ! - ML - global conservation diagnostics
!!st
!!gm : This does not only concern the dynamics ==>>> add a new title
!!gm2: why ouput restart before AGRIF update?
!!
!!jc: That would be better, but see comment above
!!
!!====>>>> to be modified for RK3
IF( lrst_oce ) CALL rst_write ( kstp, Nbb, Nnn, Naa ) ! write output ocean restart file
IF( ln_sto_eos ) CALL sto_rst_write( kstp ) ! write restart file for stochastic parameters
#if defined key_agrif
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! AGRIF recursive integration
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
Kbb_a = Nbb; Kmm_a = Nbb; Krhs_a = Nrhs ! agrif_oce module copies of time level indices
CALL Agrif_Integrate_ChildGrids( stp_RK3 ) ! allows to finish all the Child Grids before updating
#endif
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Control
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
CALL stp_ctl ( kstp, Nbb )
#if defined key_agrif
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! AGRIF update
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( Agrif_NbStepint() == 0 .AND. nstop == 0 ) &
& CALL Agrif_update_all( ) ! Update all components
#endif
IF( ln_diaobs .AND. nstop == 0 ) &
& CALL dia_obs( kstp, Nnn ) ! obs-minus-model (assimilation) diags (after dynamics update)
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! File manipulation at the end of the first time step
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( kstp == nit000 ) THEN ! 1st time step only
CALL iom_close( numror ) ! close input ocean restart file
IF( lrxios ) CALL iom_context_finalize( cr_ocerst_cxt )
IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce
IF(lwm .AND. numoni /= -1 ) CALL FLUSH ( numoni ) ! flush output namelist ice (if exist)
ENDIF
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Coupled mode
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( lk_oasis .AND. nstop == 0 ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges
!
#if defined key_xios
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Finalize contextes if end of simulation or error detected
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( kstp == nitend .OR. nstop > 0 ) THEN
CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
ENDIF
#endif
!
IF( ln_timing ) CALL timing_stop('stp_RK3')
!
END SUBROUTINE stp_RK3
SUBROUTINE mlf_baro_corr( Kmm, Kaa, puu, pvv )
!!----------------------------------------------------------------------
!! *** ROUTINE mlf_baro_corr ***
!!
!! ** Purpose : Finalize after horizontal velocity.
!!
!! ** Method : * Ensure after velocities transport matches time splitting
!! estimate (ln_dynspg_ts=T)
!!
!! ** Action : puu(Kmm),pvv(Kmm) updated now horizontal velocity (ln_bt_fw=F)
!! puu(Kaa),pvv(Kaa) after horizontal velocity
!!----------------------------------------------------------------------
!!
INTEGER , INTENT(in ) :: Kmm, Kaa ! before and after time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities
!
INTEGER :: jk ! dummy loop indices
REAL(wp), DIMENSION(jpi,jpj) :: zue, zve
!!----------------------------------------------------------------------
! Ensure below that barotropic velocities match time splitting estimate
! Compute actual transport and replace it with ts estimate at "after" time step
zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
DO jk = 2, jpkm1
zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
END DO
DO jk = 1, jpkm1
puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
END DO
!
!!st IF( .NOT.ln_bt_fw ) THEN
!!st ! Remove advective velocity from "now velocities"
!!st ! prior to asselin filtering
!!st ! In the forward case, this is done below after asselin filtering
!!st ! so that asselin contribution is removed at the same time
!!st DO jk = 1, jpkm1
!!st puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)
!!st pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
!!st END DO
!!st ENDIF
!
END SUBROUTINE mlf_baro_corr
SUBROUTINE rk3_dia( kstp, kswitch )
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kstp ! ocean time-step index
INTEGER, INTENT(in) :: kswitch ! on/off = 1/0
!!
LOGICAL, SAVE :: ll_trddyn, ll_trdtrc, ll_trdtra ! call trd at stage 3 only
LOGICAL, SAVE :: ll_diaptr, ll_ldfeiv_dia
!!----------------------------------------------------------------------
!
IF( kstp == nit000 ) THEN ! save diagnotic logical
ll_trdtra = l_trdtra
ll_trdtrc = l_trdtrc
ll_trddyn = l_trddyn
ll_diaptr = l_diaptr
ll_ldfeiv_dia = l_ldfeiv_dia
ENDIF
!
SELECT CASE( kswitch )
CASE ( 1 ) ! diagnostic activated (on)
l_trdtra = ll_trdtra
l_trdtrc = ll_trdtrc
l_trddyn = ll_trddyn
l_diaptr = ll_diaptr
l_ldfeiv_dia = ll_ldfeiv_dia
CASE ( 0 ) ! diagnostic desactivated (off)
l_trdtra = .FALSE.
l_trdtrc = .FALSE.
l_trddyn = .FALSE.
l_diaptr = .FALSE.
l_ldfeiv_dia = .FALSE.
END SELECT
!
END SUBROUTINE rk3_dia
#else
!!----------------------------------------------------------------------
!! default option EMPTY MODULE qco not activated
!!----------------------------------------------------------------------
#endif
!!======================================================================
END MODULE stprk3
src/OCE/stprk3_stg.F90 0 → 100644
View file @ d1c6ae27
MODULE stprk3_stg
!!======================================================================
!! *** MODULE stprk3_stg ***
!! Time-stepping : manager of the ocean, tracer and ice time stepping
!! using a 3rd order Runge-Kutta with fixed or quasi-eulerian coordinate
!!======================================================================
!! History : 4.5 ! 2021-01 (S. Techene, G. Madec, N. Ducousso, F. Lemarie) Original code
!! NEMO
!!----------------------------------------------------------------------
#if defined key_qco || defined key_linssh
!!----------------------------------------------------------------------
!! 'key_qco' Quasi-Eulerian vertical coordinate
!! OR
!! 'key_linssh Fixed in time vertical coordinate
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! stp_RK3_stg : NEMO 3rd order Runge-Kutta stage with qco or linssh
!!----------------------------------------------------------------------
USE step_oce ! time stepping used modules
USE domqco ! quasi-eulerian coordinate (dom_qco_r3c routine)
USE bdydyn ! ocean open boundary conditions (define bdy_dyn)
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
# if defined key_top
USE trc ! ocean passive tracers variables
USE trcadv ! passive tracers advection (trc_adv routine)
USE trcsms ! passive tracers source and sink
USE trctrp ! passive tracers transport
USE trcsbc ! passive tracers surface boundary condition !!st WARNING USELESS TO BE REMOVED
USE trcbdy ! passive tracers transport open boundary
USE trcstp_rk3
# endif
# if defined key_agrif
USE agrif_oce_interp
# endif
!
USE prtctl ! print control
IMPLICIT NONE
PRIVATE
PUBLIC stp_RK3_stg ! called by nemogcm.F90
REAL(wp) :: r1_2 = 1._wp / 2._wp
REAL(wp) :: r1_3 = 1._wp / 3._wp
REAL(wp) :: r2_3 = 2._wp / 3._wp
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssha ! sea-surface height at N+1
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ua_b, va_b ! barotropic velocity at N+1
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3ta, r3ua, r3va ! ssh/h_0 ratio at t,u,v-column at N+1
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3fb, r3fa ! ssh/h_0 ratio at f-column at N and N+1
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE stp_RK3_stg( kstg, kstp, Kbb, Kmm, Krhs, Kaa )
!!----------------------------------------------------------------------
!! *** ROUTINE stp_RK3_stg ***
!!
!! ** Purpose : - stage of RK3 time stepping of OCE and TOP
!!
!! ** Method : input: computed in dynspg_ts
!! ssh shea surface height at N+1 (oce.F90)
!! (uu_b,vv_b) barotropic velocity at N, N+1 (oce.F90)
!! (un_adv,vn_adv) barotropic transport from N to N+1 (dynspg_ts.F90)
!! ,
!! -1- set ssh(Naa) (Naa=N+1/3, N+1/2, or N)
!! -2- set the advective velocity (zadU,zaV)
!! -4- Compute the after (Naa) T-S
!! -5- Update now
!! -6- Update the horizontal velocity
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kstg ! RK3 stage
INTEGER, INTENT(in) :: kstp, Kbb, Kmm, Krhs, Kaa ! ocean time-step and time-level indices
!
INTEGER :: ji, jj, jk, jn, jtile ! dummy loop indices
REAL(wp) :: ze3Tb, ze3Sb, z1_e3t ! local scalars
REAL(wp) :: ze3Tr, ze3Sr ! - -
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zaU, zaV ! advective horizontal velocity
REAL(wp), DIMENSION(jpi,jpj) :: zub, zvb ! advective transport
!! ---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('stp_RK3_stg')
!
IF( kstp == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'stp_RK3_stg : Runge Kutta 3rd order at stage ', kstg
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
!
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! ssh, uu_b, vv_b, and ssh/h0 at Kaa
! 3D advective velocity at Kmm
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
SELECT CASE( kstg )
! !---------------!
CASE ( 1 ) !== Stage 1 ==! Kbb = Kmm = N ; Kaa = N+1/3
! !---------------!
!
ALLOCATE( ssha(jpi,jpj) , ua_b(jpi,jpj) , va_b(jpi,jpj) )
!
rDt = r1_3 * rn_Dt ! set time-step : rn_Dt/3
r1_Dt = 1._wp / rDt
!
ssha(:,:) = ssh (:,:,Kaa) ! save ssh, uu_b, vv_b at N+1 (computed in dynspg_ts)
ua_b(:,:) = uu_b(:,:,Kaa)
va_b(:,:) = vv_b(:,:,Kaa)
! ! interpolated ssh and (uu_b,vv_b) at Kaa (N+1/3)
ssh (:,:,Kaa) = r2_3 * ssh (:,:,Kbb) + r1_3 * ssha(:,:)
uu_b(:,:,Kaa) = r2_3 * uu_b(:,:,Kbb) + r1_3 * ua_b(:,:)
vv_b(:,:,Kaa) = r2_3 * vv_b(:,:,Kbb) + r1_3 * va_b(:,:)
!
!
! !== ssh/h0 ratio at Kaa ==!
!
IF( .NOT.lk_linssh ) THEN ! "after" ssh/h_0 ratio at t,u,v-column computed at N+1 stored in r3.a
!
ALLOCATE( r3ta(jpi,jpj) , r3ua(jpi,jpj) , r3va(jpi,jpj) , r3fa(jpi,jpj) , r3fb(jpi,jpj) )
!
r3fb(:,:) = r3f(:,:)
CALL dom_qco_r3c_RK3( ssha, r3ta, r3ua, r3va, r3fa )
!
CALL lbc_lnk( 'stprk3_stg', r3ua, 'U', 1._wp, r3va, 'V', 1._wp, r3fa, 'F', 1._wp )
! !
r3t(:,:,Kaa) = r2_3 * r3t(:,:,Kbb) + r1_3 * r3ta(:,:) ! at N+1/3 (Kaa)
r3u(:,:,Kaa) = r2_3 * r3u(:,:,Kbb) + r1_3 * r3ua(:,:)
r3v(:,:,Kaa) = r2_3 * r3v(:,:,Kbb) + r1_3 * r3va(:,:)
! r3f already properly set up ! at N (Kmm)
ENDIF
!
! !---------------!
CASE ( 2 ) !== Stage 2 ==! Kbb = N ; Kmm = N+1/3 ; Kaa = N+1/2
! !---------------!
!
rDt = r1_2 * rn_Dt ! set time-step : rn_Dt/2
r1_Dt = 1._wp / rDt
!
! ! set ssh and (uu_b,vv_b) at N+1/2 (Kaa)
ssh (:,:,Kaa) = r1_2 * ( ssh (:,:,Kbb) + ssha(:,:) )
uu_b(:,:,Kaa) = r1_2 * ( uu_b(:,:,Kbb) + ua_b(:,:) )
vv_b(:,:,Kaa) = r1_2 * ( vv_b(:,:,Kbb) + va_b(:,:) )
!
IF( .NOT.lk_linssh ) THEN
r3t(:,:,Kaa) = r1_2 * ( r3t(:,:,Kbb) + r3ta(:,:) ) ! at N+1/2 (Kaa)
r3u(:,:,Kaa) = r1_2 * ( r3u(:,:,Kbb) + r3ua(:,:) )
r3v(:,:,Kaa) = r1_2 * ( r3v(:,:,Kbb) + r3va(:,:) )
r3f(:,:) = r2_3 * r3fb(:,:) + r1_3 * r3fa(:,:) ! at N+1/3 (Kmm)
ENDIF
!
! !---------------!
CASE ( 3 ) !== Stage 3 ==! Kbb = N ; Kmm = N+1/2 ; Kaa = N+1
! !---------------!
!
rDt = rn_Dt ! set time-step : rn_Dt
r1_Dt = 1._wp / rDt
!
ssh (:,:,Kaa) = ssha(:,:) ! recover ssh and (uu_b,vv_b) at N + 1
uu_b(:,:,Kaa) = ua_b(:,:)
vv_b(:,:,Kaa) = va_b(:,:)
!
DEALLOCATE( ssha , ua_b , va_b )
!
IF( .NOT.lk_linssh ) THEN
r3t(:,:,Kaa) = r3ta(:,:) ! at N+1 (Kaa)
r3u(:,:,Kaa) = r3ua(:,:)
r3v(:,:,Kaa) = r3va(:,:)
r3f(:,: ) = r1_2 * ( r3fb(:,:) + r3fa(:,:) ) ! at N+1/2 (Kmm)
DEALLOCATE( r3ta, r3ua, r3va, r3fb ) ! deallocate all r3. except r3fa which will be
! ! saved in r3f at the end of the time integration and then deallocated
!
ENDIF
!
END SELECT
!
! !== advective velocity at Kmm ==!
!
! !- horizontal components -! (zaU,zaV)
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
zub(ji,jj) = un_adv(ji,jj)*r1_hu(ji,jj,Kmm) - uu_b(ji,jj,Kmm) ! barotropic velocity correction
zvb(ji,jj) = vn_adv(ji,jj)*r1_hv(ji,jj,Kmm) - vv_b(ji,jj,Kmm)
END_2D
DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) ! horizontal advective velocity
zaU(ji,jj,jk) = uu(ji,jj,jk,Kmm) + zub(ji,jj)*umask(ji,jj,jk)
zaV(ji,jj,jk) = vv(ji,jj,jk,Kmm) + zvb(ji,jj)*vmask(ji,jj,jk)
END_3D
! !- vertical components -! ww
!
IF( ln_dynadv_vec ) THEN ! ww cross-level velocity consistent with uu/vv at Kmm
CALL wzv ( kstp, Kbb, Kmm, Kaa, uu(:,:,:,Kmm), vv(:,:,:,Kmm), ww )
ELSE ! ww cross-level velocity consistent with zaU/zaV
CALL wzv ( kstp, Kbb, Kmm, Kaa, zaU, zaV, ww )
ENDIF
!!st IF( ln_zad_Aimp ) CALL wAimp( kstp, ww, wi ) ! Adaptive-implicit vertical advection partitioning
! !==>>> implicite for stages 1 & 2 ????
!
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RHS of ocean dynamics : ADV + VOR/COR + HPG (+ ASM ) <<<=== Question: Stokes drift ?
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
uu(:,:,:,Krhs) = 0._wp ! set dynamics trends to zero
vv(:,:,:,Krhs) = 0._wp
!
!===>>>>>> Modify dyn_adv_... dyn_keg routines so that Krhs to zero useless
! ! advection (VF or FF) ==> RHS
IF( ln_dynadv_vec ) THEN ! uu and vv used for momentum advection
CALL dyn_adv( kstp, Kbb, Kmm , uu, vv, Krhs)
ELSE ! advective velocity used for momentum advection
CALL dyn_adv( kstp, Kbb, Kmm , uu, vv, Krhs, zaU, zaV, ww )
ENDIF
! ! Coriolis / vorticity ==> RHS
CALL dyn_vor( kstp, Kmm , uu, vv, Krhs )
!
!!gm à appeler que pour ln_zad_Aimp=T et en ne faisant que wi par zdf
!! ! ZAD (implicit part) ==> RHS
!! CALL dyn_zdf ( kstp, Kbb, Kmm, Krhs, uu, vv, Kaa )
!===>>>>>> Modify dyn_hpg & dyn_hpg_... routines : rhd computed in dyn_hpg and pass in argument to dyn_hpg_...
CALL eos ( ts, Kmm, rhd ) ! Kmm in situ density anomaly for hpg computation
!!gm end
CALL dyn_hpg( kstp, Kmm , uu, vv, Krhs )
!
!!gm ===>>>>>> Probably useless since uu_b(Kaa) will be imposed at the end of stage 1 and 2
! but may be necessary in stage 3 due to implicite in dynzdf.
! except if my idea for the matrice construction is OK !
! ! ! grad_h of ps ==> RHS
! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
! uu(ji,jj,jk,Krhs) = uu(ji,jj,jk,Krhs) - grav * ( ssh(ji+1,jj ,Kmm) - ssh(ji,jj,Kmm) )
! vv(ji,jj,jk,Krhs) = vv(ji,jj,jk,Krhs) - grav * ( ssh(ji ,jj+1,Kmm) - ssh(ji,jj,Kmm) )
! END_3D
!!gm
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RHS of tracers : ADV only using (zaU,zaV,ww)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
! ! Advective velocity needed for tracers advection - already computed if ln_dynadv_vec=F
IF( ln_dynadv_vec ) CALL wzv ( kstp, Kbb, Kmm, Kaa, zaU, zaV, ww )
!
# if defined key_top
! !== Passive Tracer ==!
!
SELECT CASE( kstg )
! !-------------------!
CASE ( 1 , 2 ) !== Stage 1 & 2 ==! stg1: Kbb = N ; Kaa = N+1/3
! !-------------------! stg2: Kbb = N ; Kmm = N+1/3 ; Kaa = N+1/2
!
IF( kstg == 1 ) THEN
CALL trc_stp_start( kstp, Kbb, Kmm, Krhs, Kaa )
ENDIF
!
IF(.NOT. ln_trcadv_mus ) THEN
!
DO jn = 1, jptra
tr(:,:,:,jn,Krhs) = 0._wp ! set tracer trends to zero !!st ::: required because of tra_adv new loops
END DO
! !== advection of passive tracers ==!
rDt_trc = rDt
!
CALL trc_sbc_RK3( kstp, Kmm, tr, Krhs, kstg ) ! surface boundary condition
!
CALL trc_adv ( kstp, Kbb, Kmm, tr, Krhs, zaU, zaV, ww ) ! horizontal & vertical advection
!
! !== time integration ==! ∆t = rn_Dt/3 (stg1) or rn_Dt/2 (stg2)
DO jn = 1, jptra
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
ze3Tb = e3t(ji,jj,jk,Kbb) * tr(ji,jj,jk,jn,Kbb )
ze3Tr = e3t(ji,jj,jk,Kmm) * tr(ji,jj,jk,jn,Krhs)
z1_e3t= 1._wp / e3t(ji,jj,jk, Kaa)
tr(ji,jj,jk,jn,Kaa) = ( ze3Tb + rDt * ze3Tr*tmask(ji,jj,jk) ) * z1_e3t
END_3D
END DO
!
!!st need a lnc lkn at stage 1 & 2 otherwise tr@Kmm will not be usable in trc_adv
CALL lbc_lnk( 'stprk3_stg', tr(:,:,:,:,Kaa), 'T', 1._wp )
ENDIF
! !---------------!
CASE ( 3 ) !== Stage 3 ==! add all RHS terms but advection (=> Kbb only)
! !---------------!
!
DO jn = 1, jptra
tr(:,:,:,jn,Krhs) = 0._wp
END DO
! !== advection of passive tracers ==!
rDt_trc = rDt
!
CALL trc_sbc_RK3( kstp, Kmm, tr, Krhs, kstg ) ! surface boundary condition
!
CALL trc_adv ( kstp, Kbb, Kmm, tr, Krhs, zaU, zaV, ww ) ! horizontal & vertical advection
!
CALL trc_sms ( kstp, Kbb, Kbb, Krhs ) ! tracers: sinks and sources
CALL trc_trp ( kstp, Kbb, Kmm, Krhs, Kaa ) ! transport of passive tracers (without advection)
!
!
CALL trc_stp_end( kstp, Kbb, Kmm, Kaa )
!
END SELECT
# endif
! !== T-S Tracers ==!
!===>>>>>> Modify tra_adv_... routines so that Krhs to zero useless
DO jn = 1, jpts
ts(:,:,:,jn,Krhs) = 0._wp ! set tracer trends to zero (:,:,:) needed otherwise it does not work (?)
END DO
CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0 ) ! now in potential density for tra_mle computation
!===>>> CAUTION here may be without GM velocity but stokes drift required ! 0 barotropic divergence for GM != 0 barotropic divergence for SD
!!st consistence 2D / 3D - flux de masse
CALL tra_adv( kstp, Kbb, Kmm, ts, Krhs, zaU, zaV, ww ) ! hor. + vert. advection ==> RHS
!===>>>>>> stg1&2: Verify the necessity of these trends (we may need it as there are in the RHS of dynspg_ts ?)
!!gm ====>>>> needed for heat and salt fluxes associated with mass/volume flux
CALL tra_sbc_RK3( kstp, Kmm, ts, Krhs, kstg ) ! surface boundary condition
IF( ln_isf ) CALL tra_isf ( kstp, Kmm, ts, Krhs ) ! ice shelf heat flux
IF( ln_traqsr ) CALL tra_qsr ( kstp, Kmm, ts, Krhs ) ! penetrative solar radiation qsr
!!gm
!
!!gm ===>>>>>> Verify the necessity of these trends at stages 1 and 2
! (we may need it as they are in the RHS of dynspg_ts ?)
! IF( lk_asminc .AND. ln_asmiau ) THEN ! apply assimilation increment
! IF( ln_dyninc ) CALL dyn_asm_inc( kstp, Kbb, Kmm, uu, vv, Krhs ) ! dynamics ==> RHS
! IF( ln_trainc ) CALL tra_asm_inc( kstp, Kbb, Kmm, ts , Krhs ) ! tracers ==> RHS
! ENDIF
!!gm end Verif
!
SELECT CASE( kstg )
! !-------------------!
CASE ( 1 , 2 ) !== Stage 1 & 2 ==! stg1: Kbb = N ; Kaa = N+1/3
! !-------------------! stg2: Kbb = N ; Kmm = N+1/3 ; Kaa = N+1/2
!
! !== time integration ==! ∆t = rn_Dt/3 (stg1) or rn_Dt/2 (stg2)
IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
uu(ji,jj,jk,Kaa) = ( uu(ji,jj,jk,Kbb) + rDt * uu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
vv(ji,jj,jk,Kaa) = ( vv(ji,jj,jk,Kbb) + rDt * vv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
END_3D
ELSE ! applied on thickness weighted velocity
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
uu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb ) &
& + rDt * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Krhs) ) &
& / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
vv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb ) &
& + rDt * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Krhs) ) &
& / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
END_3D
ENDIF
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
ze3Tb = e3t(ji,jj,jk,Kbb) * ts(ji,jj,jk,jp_tem,Kbb )
ze3Sb = e3t(ji,jj,jk,Kbb) * ts(ji,jj,jk,jp_sal,Kbb )
ze3Tr = e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Krhs)
ze3Sr = e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Krhs)
z1_e3t= 1._wp / e3t(ji,jj,jk, Kaa)
ts(ji,jj,jk,jp_tem,Kaa) = ( ze3Tb + rDt * ze3Tr*tmask(ji,jj,jk) ) * z1_e3t
ts(ji,jj,jk,jp_sal,Kaa) = ( ze3Sb + rDt * ze3Sr*tmask(ji,jj,jk) ) * z1_e3t
END_3D
!
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=uu(:,:,:,Kaa), clinfo1='stp stg - Ua: ', mask1=umask, &
& tab3d_2=vv(:,:,:,Kaa), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Kaa), clinfo1='stp stg - Ta: ', mask1=tmask, &
& tab3d_2=ts(:,:,:,jp_sal,Kaa), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
!
! !---------------!
CASE ( 3 ) !== Stage 3 ==! add all remaining RHS terms
! !---------------!
!
! !== complete the momentum RHS ==! except ZDF (implicit)
! ! lateral mixing ==> RHS
CALL dyn_ldf( kstp, Kbb, Kmm, uu, vv, Krhs )
! ! OSMOSIS non-local velocity fluxes ==> RHS
IF( ln_zdfosm ) CALL dyn_osm( kstp, Kmm, uu, vv, Krhs )
!
IF( ln_bdy ) THEN ! bdy damping trends ==> RHS
CALL bdy_dyn3d_dmp ( kstp, Kbb, uu, vv, Krhs )
CALL bdy_tra_dmp ( kstp, Kbb, ts , Krhs )
ENDIF
# if defined key_agrif
IF(.NOT. Agrif_Root() ) THEN ! AGRIF: sponge ==> momentum and tracer RHS
CALL Agrif_Sponge_dyn
CALL Agrif_Sponge_tra
ENDIF
# endif
! !== complete the tracers RHS ==! except ZDF (implicit)
! !* T-S Tracer *!
!
CALL tra_ldf( kstp, Kbb, Kmm, ts, Krhs ) ! lateral mixing
IF( ln_trabbc ) CALL tra_bbc( kstp, Kmm, ts, Krhs ) ! bottom heat flux
IF( ln_trabbl ) CALL tra_bbl( kstp, Kbb, Kmm, ts, Krhs ) ! advective (and/or diffusive) bottom boundary layer scheme
IF( ln_tradmp ) CALL tra_dmp( kstp, Kbb, Kmm, ts, Krhs ) ! internal damping trends
IF( ln_zdfmfc ) CALL tra_mfc( kstp, Kbb, ts, Krhs ) ! Mass Flux Convection
IF( ln_zdfosm ) THEN
CALL tra_osm( kstp, Kmm, ts, Krhs ) ! OSMOSIS non-local tracer fluxes ==> RHS
IF( lrst_oce ) CALL osm_rst( kstp, Kmm, 'WRITE' ) ! write OSMOSIS outputs + ww (so must do here) to restarts
ENDIF
!
! !== DYN & TRA time integration + ZDF ==! ∆t = rDt
!
CALL dyn_zdf( kstp, Kbb, Kmm, Krhs, uu, vv, Kaa ) ! vertical diffusion and time integration
CALL tra_zdf( kstp, Kbb, Kmm, Krhs, ts , Kaa ) ! vertical mixing and after tracer fields
IF( ln_zdfnpc ) CALL tra_npc( kstp, Kmm, Krhs, ts , Kaa ) ! update after fields by non-penetrative convection
!
IF( .NOT.lk_linssh ) THEN
r3f(:,:) = r3fa(:,:) ! save r3fa in r3f before deallocation
DEALLOCATE( r3fa ) ! (r3f = r3f(Kbb) of the next time step)
ENDIF
!
END SELECT
! !== correction of the barotropic (all stages) ==! at Kaa = N+1/3, N+1/2 or N+1
! ! barotropic velocity correction
zub(A2D(0)) = uu_b(A2D(0),Kaa) - SUM( e3u_0(A2D(0),:)*uu(A2D(0),:,Kaa), 3 ) * r1_hu_0(A2D(0))
zvb(A2D(0)) = vv_b(A2D(0),Kaa) - SUM( e3v_0(A2D(0),:)*vv(A2D(0),:,Kaa), 3 ) * r1_hv_0(A2D(0))
!
DO jk = 1, jpkm1 ! corrected horizontal velocity
uu(A2D(0),jk,Kaa) = uu(A2D(0),jk,Kaa) + zub(A2D(0))*umask(A2D(0),jk)
vv(A2D(0),jk,Kaa) = vv(A2D(0),jk,Kaa) + zvb(A2D(0))*vmask(A2D(0),jk)
END DO
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Set boundary conditions
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
# if defined key_agrif
CALL Agrif_tra( kstp, kstg ) !* AGRIF zoom boundaries
CALL Agrif_dyn( kstp, kstg )
# endif
! !* local domain boundaries (T-point, unchanged sign)
CALL lbc_lnk( 'stp_RK3_stg', uu(:,:,:, Kaa), 'U', -1., vv(:,:,: ,Kaa), 'V', -1. &
& , ts(:,:,:,jp_tem,Kaa), 'T', 1., ts(:,:,:,jp_sal,Kaa), 'T', 1. )
!
IF( nn_hls == 2 .AND. l_zdfsh2 ) CALL lbc_lnk( 'stp_RK3_stg', avm_k, 'W', 1.0_wp ) ! lbc_lnk needed for zdf_sh2 when using nn_hls = 2, moved here to allow tiling in zdf_phy
!
! !* BDY open boundaries
IF( ln_bdy ) THEN
CALL bdy_tra( kstp, Kbb, ts, Kaa )
IF( ln_dynspg_exp ) CALL bdy_dyn( kstp, Kbb, uu, vv, Kaa )
IF( ln_dynspg_ts ) CALL bdy_dyn( kstp, Kbb, uu, vv, Kaa, dyn3d_only=.true. )
ENDIF
!
IF( ln_timing ) CALL timing_stop('stp_RK3_stg')
!
END SUBROUTINE stp_RK3_stg
#else
!!----------------------------------------------------------------------
!! default option EMPTY MODULE qco not activated
!!----------------------------------------------------------------------
#endif
!!======================================================================
END MODULE stprk3_stg
src/OFF/dtadyn.F90
View file @ d1c6ae27
......@@ -101,7 +101,7 @@ MODULE dtadyn
!!----------------------------------------------------------------------
!! NEMO/OFF 4.0 , NEMO Consortium (2018)
!! $Id: dtadyn.F90 15090 2021-07-06 14:25:18Z cetlod $
!! $Id: dtadyn.F90 15532 2021-11-24 11:47:32Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -175,7 +175,7 @@ CONTAINS
ENDIF
ENDIF
!
CALL eos ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
CALL eos ( ts(:,:,:,:,Kmm), rhd, gdept_0(:,:,:) ) ! In any case, we need rhd
CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm ) ! now local thermal/haline expension ratio at T-points
CALL bn2 ( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm ) ! before Brunt-Vaisala frequency need for zdfmxl
......@@ -192,7 +192,7 @@ CONTAINS
ENDIF
!
!
CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
CALL eos( ts(:,:,:,:,Kmm), rhd, gdept_0(:,:,:) ) ! In any case, we need rhd
!
IF(sn_cfctl%l_prtctl) THEN ! print control
CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn - : ', mask1=tmask )
......@@ -669,7 +669,7 @@ CONTAINS
!!---------------------------------------------------------------------
!
IF( l_ldfslp .AND. .NOT.ln_c1d ) THEN ! Computes slopes (here avt is used as workspace)
CALL eos ( pts, rhd, rhop, gdept_0(:,:,:) )
CALL eos ( pts, rhd, gdept_0(:,:,:) )
CALL eos_rab( pts, rab_n, Kmm ) ! now local thermal/haline expension ratio at T-points
CALL bn2 ( pts, rab_n, rn2, Kmm ) ! now Brunt-Vaisala
......@@ -727,7 +727,7 @@ CONTAINS
ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature
ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity
!
CALL eos ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
CALL eos ( ts(:,:,:,:,Kmm), rhd, gdept_0(:,:,:) ) ! In any case, we need rhd
IF(sn_cfctl%l_prtctl) THEN ! print control
CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn - : ', mask1=tmask )
......
src/OFF/nemogcm.F90
View file @ d1c6ae27
......@@ -82,7 +82,7 @@ MODULE nemogcm
!!----------------------------------------------------------------------
!! NEMO/OFF 4.0 , NEMO Consortium (2018)
!! $Id: nemogcm.F90 15446 2021-10-26 14:34:38Z cetlod $
!! $Id: nemogcm.F90 14255 2021-01-04 11:35:00Z cetlod $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -517,7 +517,6 @@ CONTAINS
ts (:,:,:,:,Kmm) = 0._wp ! !
!
rhd (:,:,:) = 0.e0
rhop (:,:,:) = 0.e0
rn2 (:,:,:) = 0.e0
!
END SUBROUTINE istate_init
......
src/TOP/AGE/trcsms_age.F90
View file @ d1c6ae27
......@@ -31,7 +31,7 @@ MODULE trcsms_age
!!----------------------------------------------------------------------
!! NEMO/TOP 4.0 , NEMO Consortium (2018)
!! $Id: trcsms_age.F90 14173 2020-12-15 12:44:07Z cetlod $
!! $Id: trcsms_age.F90 15193 2021-08-13 13:18:24Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -51,14 +51,17 @@ CONTAINS
!
IF( ln_timing ) CALL timing_start('trc_sms_age')
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' trc_sms_age: AGE model'
IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
IF( kt == nittrc000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' trc_sms_age: AGE model'
IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
ENDIF
#if ! defined key_RK3
IF( l_1st_euler .OR. ln_top_euler ) THEN
tr(:,:,:,jp_age,Kbb) = tr(:,:,:,jp_age,Kmm)
ENDIF
#endif
DO jk = 1, nla_age
tr(:,:,jk,jp_age,Krhs) = rn_age_kill_rate * tr(:,:,jk,jp_age,Kbb)
......@@ -71,7 +74,7 @@ CONTAINS
tr(:,:,jk,jp_age,Krhs) = tmask(:,:,jk) * rryear
END DO
!
IF( l_trdtrc ) CALL trd_trc( tr(:,:,:,jp_age,Krhs), jn, jptra_sms, kt, Kmm ) ! save trends
IF( l_trdtrc ) CALL trd_trc( tr(:,:,:,jp_age,Krhs), jn, jptra_sms, kt, Kmm ) ! save trends
!
IF( ln_timing ) CALL timing_stop('trc_sms_age')
!
......
src/TOP/PISCES/P4Z/p4zche.F90
View file @ d1c6ae27
......@@ -11,6 +11,7 @@ MODULE p4zche
!! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90
!! ! 2011-02 (J. Simeon, J.Orr ) update O2 solubility constants
!! 3.6 ! 2016-03 (O. Aumont) Change chemistry to MOCSY standards
!! 4.2 ! 2020 (J. ORR ) rhop is replaced by "in situ density" rhd
!!----------------------------------------------------------------------
!! p4z_che : Sea water chemistry computed following OCMIP protocol
!!----------------------------------------------------------------------
......@@ -134,7 +135,7 @@ MODULE p4zche
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/TOP 4.0 , NEMO Consortium (2018)
!! $Id: p4zche.F90 15459 2021-10-29 08:19:18Z cetlod $
!! $Id: p4zche.F90 15532 2021-11-24 11:47:32Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -194,9 +195,21 @@ CONTAINS
zsal = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
! ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
! ! AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
! J. ORR: The previous code has been modified. It computed CO2 solubility in mol/(kg*atm), then converted that to mol/(L*atm).
! But Weiss (1974) provides sets of coefficients for each of those 2 units.
! Thus I have changed the code to use the coefficients for mol*L/atm.
! Hence I've eliminated using the conversion (which used the variable rhop)
! OLD - Coefficients for CO2 soulbility in mol/(kg*atm) (Weiss,1974, Table 1, column 2)
!zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel &
!& + 0.0047036e-4*ztkel**2)
! NEW - Coefficients for CO2 soulbility in mol/(L*atm) (Weiss, 1974, Table 1, column 1)
zcek1 = 9050.69/ztkel - 58.0931 + 22.2940 * LOG(zt) + zsal*(0.027766 - 0.00025888*ztkel &
& + 0.0050578e-4*ztkel**2)
chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 ! mol/(L atm)
& + 0.0050578e-4*ztkel**2)
!
! OLD: chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm)
! The units indicated in the above line are wrong. They are actually "mol/(L*uatm)"
! NEW:
chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 ! mol/(L * uatm)
chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3
chemc(ji,jj,3) = 57.7 - 0.118*ztkel
END_2D
......@@ -444,15 +457,16 @@ CONTAINS
INTEGER :: ji, jj, jk
REAL(wp) :: zca1, zba1
REAL(wp) :: zd, zsqrtd, zhmin
REAL(wp) :: za2, za1, za0
REAL(wp) :: za2, za1, za0, zrhd
REAL(wp) :: p_dictot, p_bortot, p_alkcb
!!---------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('ahini_for_at')
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
p_alkcb = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
p_dictot = tr(ji,jj,jk,jpdic,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
p_alkcb = tr(ji,jj,jk,jptal,Kbb) * zrhd
p_dictot = tr(ji,jj,jk,jpdic,Kbb) * zrhd
p_bortot = borat(ji,jj,jk)
IF (p_alkcb <= 0.) THEN
p_hini(ji,jj,jk) = 1.e-3
......@@ -501,11 +515,16 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup
INTEGER, INTENT(in) :: Kbb ! time level indices
INTEGER :: ji, jj, jk
REAL(wp) :: zrhd
p_alknw_inf(:,:,:) = -tr(:,:,:,jppo4,Kbb) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:) &
& - fluorid(:,:,:)
p_alknw_sup(:,:,:) = (2. * tr(:,:,:,jpdic,Kbb) + 2. * tr(:,:,:,jppo4,Kbb) + tr(:,:,:,jpsil,Kbb) ) &
& * 1000. / (rhop(:,:,:) + rtrn) + borat(:,:,:)
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
p_alknw_inf(ji,jj,jk) = -tr(ji,jj,jk,jppo4,Kbb) * zrhd - sulfat(ji,jj,jk) &
& - fluorid(ji,jj,jk)
p_alknw_sup(ji,jj,jk) = (2. * tr(ji,jj,jk,jpdic,Kbb) + 2. * tr(ji,jj,jk,jppo4,Kbb) + tr(ji,jj,jk,jpsil,Kbb) ) &
& * zrhd + borat(ji,jj,jk)
END_3D
END SUBROUTINE anw_infsup
......@@ -535,7 +554,7 @@ CONTAINS
REAL(wp) :: znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4
REAL(wp) :: znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu
REAL(wp) :: zalk_wat, zdalk_wat
REAL(wp) :: zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit
REAL(wp) :: zrhd, p_alktot, zdic, zbot, zpt, zst, zft, zsit
LOGICAL :: l_exitnow
REAL(wp), PARAMETER :: pz_exp_threshold = 1.0
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin
......@@ -550,7 +569,8 @@ CONTAINS
! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
IF (rmask(ji,jj,jk) == 1.) THEN
p_alktot = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
p_alktot = tr(ji,jj,jk,jptal,Kbb) * zrhd
aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk)
zh_ini = p_hini(ji,jj,jk)
......@@ -579,12 +599,12 @@ CONTAINS
DO jn = 1, jp_maxniter_atgen
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
IF (rmask(ji,jj,jk) == 1.) THEN
zfact = rhop(ji,jj,jk) / 1000. + rtrn
p_alktot = tr(ji,jj,jk,jptal,Kbb) / zfact
zdic = tr(ji,jj,jk,jpdic,Kbb) / zfact
zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
p_alktot = tr(ji,jj,jk,jptal,Kbb) * zrhd
zdic = tr(ji,jj,jk,jpdic,Kbb) * zrhd
zbot = borat(ji,jj,jk)
zpt = tr(ji,jj,jk,jppo4,Kbb) / zfact * po4r
zsit = tr(ji,jj,jk,jpsil,Kbb) / zfact
zpt = tr(ji,jj,jk,jppo4,Kbb) * zrhd * po4r
zsit = tr(ji,jj,jk,jpsil,Kbb) * zrhd
zst = sulfat (ji,jj,jk)
zft = fluorid(ji,jj,jk)
aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk)
......
src/TOP/PISCES/P4Z/p4zflx.F90
View file @ d1c6ae27
......@@ -9,6 +9,7 @@ MODULE p4zflx
!! 1.0 ! 2004 (O. Aumont) modifications
!! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90
!! ! 2011-02 (J. Simeon, J. Orr) Include total atm P correction
!! 4.2 ! 2020 (J. ORR ) rhop is replaced by "in situ density" rhd
!!----------------------------------------------------------------------
!! p4z_flx : CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
!! p4z_flx_init : Read the namelist
......@@ -56,7 +57,7 @@ MODULE p4zflx
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/TOP 4.0 , NEMO Consortium (2018)
!! $Id: p4zflx.F90 15459 2021-10-29 08:19:18Z cetlod $
!! $Id: p4zflx.F90 15532 2021-11-24 11:47:32Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -77,7 +78,7 @@ CONTAINS
!
INTEGER :: ji, jj, jm, iind, iindm1
REAL(wp) :: ztc, ztc2, ztc3, ztc4, zws, zkgwan
REAL(wp) :: zfld, zflu, zfld16, zflu16, zfact
REAL(wp) :: zfld, zflu, zfld16, zflu16, zrhd
REAL(wp) :: zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff
REAL(wp) :: zph, zdic, zsch_o2, zsch_co2
REAL(wp) :: zyr_dec, zdco2dt
......@@ -111,9 +112,9 @@ CONTAINS
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! DUMMY VARIABLES FOR DIC, H+, AND BORATE
zfact = rhop(ji,jj,1) / 1000. + rtrn
zrhd = rhd(ji,jj,1) + 1._wp
zdic = tr(ji,jj,1,jpdic,Kbb)
zph = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
zph = MAX( hi(ji,jj,1), 1.e-10 ) / ( zrhd + rtrn )
! CALCULATE [H2CO3]
zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)
END_2D
......
src/TOP/PISCES/P4Z/p4zlys.F90
View file @ d1c6ae27
......@@ -11,6 +11,7 @@ MODULE p4zlys
!! ! 2011-02 (J. Simeon, J. Orr) Calcon salinity dependence
!! 3.4 ! 2011-06 (O. Aumont, C. Ethe) Improvment of calcite dissolution
!! 3.6 ! 2015-05 (O. Aumont) PISCES quota
!! 4.2 ! 2020 (J. ORR ) rhop is replaced by "in situ density" rhd
!!----------------------------------------------------------------------
!! p4z_lys : Compute the CaCO3 dissolution
!! p4z_lys_init : Read the namelist parameters
......@@ -32,13 +33,17 @@ MODULE p4zlys
REAL(wp), PUBLIC :: nca !: order of reaction for calcite dissolution
INTEGER :: rmtss ! number of seconds per month
REAL(wp) :: calcon = 1.03E-2 ! mean calcite concentration [Ca2+] in sea water [mole/kg solution]
!! * Module variables
REAL(wp) :: calcon = 1.03E-2 !: mean calcite concentration [Ca2+] in sea water [mole/kg solution]
! J. ORR: Made consistent with mocsy's choice based on literature review from Munhoven
! REAL(wp) :: calcon = 1.0287E-2 !: mean calcite concentration [Ca2+] in sea water [mole/kg solution]
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/TOP 4.0 , NEMO Consortium (2018)
!! $Id: p4zlys.F90 15287 2021-09-24 11:11:02Z cetlod $
!! $Id: p4zlys.F90 15532 2021-11-24 11:47:32Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
......@@ -58,7 +63,7 @@ CONTAINS
INTEGER, INTENT(in) :: Kbb, Krhs ! time level indices
!
INTEGER :: ji, jj, jk, jn
REAL(wp) :: zdispot, zfact, zcalcon
REAL(wp) :: zdispot, zrhd, zcalcon
REAL(wp) :: zomegaca, zexcess, zexcess0, zkd
CHARACTER (len=25) :: charout
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zco3, zcaldiss, zhinit, zhi, zco3sat
......@@ -66,7 +71,8 @@ CONTAINS
!
IF( ln_timing ) CALL timing_start('p4z_lys')
!
zhinit (:,:,:) = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn )
zhinit (:,:,:) = hi(:,:,:) / ( rhd(:,:,:) + 1._wp )
!
! -------------------------------------------
! COMPUTE [CO3--] and [H+] CONCENTRATIONS
......@@ -76,7 +82,7 @@ CONTAINS
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
zco3(ji,jj,jk) = tr(ji,jj,jk,jpdic,Kbb) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2 &
& + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn )
hi (ji,jj,jk) = zhi(ji,jj,jk) * rhop(ji,jj,jk) / 1000.
hi (ji,jj,jk) = zhi(ji,jj,jk) * ( rhd(ji,jj,jk) + 1._wp )
END_3D
! ---------------------------------------------------------
......@@ -88,11 +94,11 @@ CONTAINS
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
! DEVIATION OF [CO3--] FROM SATURATION VALUE
! Salinity dependance in zomegaca and divide by rhop/1000 to have good units
! Salinity dependance in zomegaca and divide by rhd to have good units
zcalcon = calcon * ( salinprac(ji,jj,jk) / 35._wp )
zfact = rhop(ji,jj,jk) / 1000._wp
zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zfact + rtrn )
zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zfact / ( zcalcon + rtrn )
zrhd = rhd(ji,jj,jk) + 1._wp
zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zrhd + rtrn )
zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zrhd / ( zcalcon + rtrn )
! SET DEGREE OF UNDER-/SUPERSATURATION
excess(ji,jj,jk) = 1._wp - zomegaca
......
src/TOP/TRP/trcadv.F90
View file @ d1c6ae27
......@@ -7,6 +7,7 @@ MODULE trcadv
!! 3.0 ! 2010-06 (C. Ethe) Adapted to passive tracers
!! 3.7 ! 2014-05 (G. Madec, C. Ethe) Add 2nd/4th order cases for CEN and FCT schemes
!! 4.0 ! 2017-09 (G. Madec) remove vertical time-splitting option
!! 4.5 ! 2021-08 (G. Madec, S. Techene) add advective velocities as optional arguments
!!----------------------------------------------------------------------
#if defined key_top
!!----------------------------------------------------------------------
......@@ -35,7 +36,7 @@ MODULE trcadv
IMPLICIT NONE
PRIVATE
PUBLIC trc_adv ! called by trctrp.F90
PUBLIC trc_adv ! called by trctrp.F90 and stprk3_stg.F90
PUBLIC trc_adv_ini ! called by trcini.F90
! !!* Namelist namtrc_adv *
......@@ -44,7 +45,7 @@ MODULE trcadv
INTEGER :: nn_cen_h, nn_cen_v ! =2/4 : horizontal and vertical choices of the order of CEN scheme
LOGICAL :: ln_trcadv_fct ! FCT scheme flag
INTEGER :: nn_fct_h, nn_fct_v ! =2/4 : horizontal and vertical choices of the order of FCT scheme
LOGICAL :: ln_trcadv_mus ! MUSCL scheme flag
LOGICAL, PUBLIC :: ln_trcadv_mus ! MUSCL scheme flag
LOGICAL :: ln_mus_ups ! use upstream scheme in vivcinity of river mouths
LOGICAL :: ln_trcadv_ubs ! UBS scheme flag
INTEGER :: nn_ubs_v ! =2/4 : vertical choice of the order of UBS scheme
......@@ -58,31 +59,35 @@ MODULE trcadv
INTEGER, PARAMETER :: np_MUS = 3 ! MUSCL scheme
INTEGER, PARAMETER :: np_UBS = 4 ! 3rd order Upstream Biased Scheme
INTEGER, PARAMETER :: np_QCK = 5 ! QUICK scheme
!! * Substitutions
# include "do_loop_substitute.h90"
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/TOP 4.0 , NEMO Consortium (2018)
!! $Id: trcadv.F90 15073 2021-07-02 14:20:14Z clem $
!! $Id: trcadv.F90 15510 2021-11-15 15:33:37Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE trc_adv( kt, Kbb, Kmm, ptr, Krhs )
SUBROUTINE trc_adv( kt, Kbb, Kmm, ptr, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE trc_adv ***
!!
!! ** Purpose : compute the ocean tracer advection trend.
!!
!! ** Method : - Update after tracers (tr(Krhs)) with the advection term following nadv
!! ** Method : - Update tr(Krhs) with the advective trend following nadv
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! ocean time-step index
INTEGER , INTENT(in) :: Kbb, Kmm, Krhs ! time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! passive tracers and RHS of tracer equation
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! time level indices
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt) , INTENT(inout) :: ptr ! passive tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk ! dummy loop index
CHARACTER (len=22) :: charout
REAL(wp), DIMENSION(:,:,:), POINTER :: zptu, zptv, zptw
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zuu, zvv, zww ! effective velocity
!!----------------------------------------------------------------------
!
......@@ -105,21 +110,32 @@ CONTAINS
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zww(ji,jj,jpk) = 0._wp
END_2D
!
IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (uu,vv,ww)
zptu => pau(:,:,:)
zptv => pav(:,:,:)
zptw => paw(:,:,:)
ELSE ! MLF: advective velocity = (uu,vv,ww)
zptu => uu(:,:,:,Kmm)
zptv => vv(:,:,:,Kmm)
zptw => ww(:,:,: )
ENDIF
!
IF( ln_wave .AND. ln_sdw ) THEN
DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) ! eulerian transport + Stokes Drift
zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * ( uu(ji,jj,jk,Kmm) + usd(ji,jj,jk) )
zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * ( vv(ji,jj,jk,Kmm) + vsd(ji,jj,jk) )
zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * ( zptu(ji,jj,jk) + usd(ji,jj,jk) )
zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * ( zptv(ji,jj,jk) + vsd(ji,jj,jk) )
END_3D
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
zww(ji,jj,jk) = e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) )
zww(ji,jj,jk) = e1e2t(ji,jj) * ( zptw(ji,jj,jk) + wsd(ji,jj,jk) )
END_3D
ELSE
DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )
zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) ! eulerian transport
zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm)
zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * zptu(ji,jj,jk) ! eulerian transport
zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * zptv(ji,jj,jk)
END_3D
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
zww(ji,jj,jk) = e1e2t(ji,jj) * ww(ji,jj,jk)
zww(ji,jj,jk) = e1e2t(ji,jj) * zptw(ji,jj,jk)
END_3D
ENDIF
!
......
src/TOP/TRP/trcatf.F90
View file @ d1c6ae27
......@@ -20,7 +20,7 @@ MODULE trcatf
!! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC
!! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename trcnxt.F90 -> trcatf.F90. Now only does time filtering.
!!----------------------------------------------------------------------
#if defined key_top
#if defined key_top && ! defined key_RK3
!!----------------------------------------------------------------------
!! 'key_top' TOP models
!!----------------------------------------------------------------------
......@@ -57,7 +57,7 @@ MODULE trcatf
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/TOP 4.0 , NEMO Consortium (2018)
!! $Id: trcatf.F90 15090 2021-07-06 14:25:18Z cetlod $
!! $Id: trcatf.F90 15373 2021-10-14 17:01:57Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -99,7 +99,7 @@ CONTAINS
ENDIF
!
#if defined key_agrif
CALL Agrif_trc ! AGRIF zoom boundaries
CALL Agrif_trc( kt ) ! AGRIF zoom boundaries
#endif
! Update after tracer on domain lateral boundaries
CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1._wp )
......
src/TOP/TRP/trcrad.F90
View file @ d1c6ae27
......@@ -34,7 +34,7 @@ MODULE trcrad
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/TOP 4.0 , NEMO Consortium (2018)
!! $Id: trcrad.F90 13324 2020-07-17 19:47:48Z acc $
!! $Id: trcrad.F90 15561 2021-11-30 13:25:02Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -148,9 +148,14 @@ CONTAINS
IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) )
zs2rdt = 1. / ( 2. * rn_Dt )
!
#if ! defined key_RK3
DO jt = 1,2 ! Loop over time indices since exactly the same fix is applied to "now" and "after" fields
IF( jt == 1 ) itime = Kbb
IF( jt == 2 ) itime = Kmm
#else
DO jt = 1,1 ! Loop over time indices since exactly the same fix is applied to "now" and "after" fields
IF( jt == 1 ) itime = Kmm
#endif
IF( PRESENT( cpreserv ) ) THEN !== total tracer concentration is preserved ==!
!
......