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 461 additions and 121 deletions
src/OCE/DOM/closea.F90
View file @ 672ff863
......@@ -258,6 +258,6 @@ CONTAINS
ALLOCATE( kmask(jpi,jpj) , STAT=ierr )
IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'alloc_csmask: failed to allocate surf array')
!
END SUBROUTINE
END SUBROUTINE alloc_csmask
END MODULE closea
src/OCE/DOM/daymod.F90
View file @ 672ff863
......@@ -68,7 +68,7 @@ CONTAINS
!! - nmonth_len, nyear_len, nmonth_beg through day_mth
!!----------------------------------------------------------------------
INTEGER :: inbday, imonday, isecrst ! local integers
REAL(wp) :: zjul ! local scalar
REAL(dp) :: zjul ! local scalar
!!----------------------------------------------------------------------
!
! max number of seconds between each restart
......@@ -93,7 +93,7 @@ CONTAINS
nminute = ( nn_time0 - nhour * 100 )
isecrst = ( nhour * NINT(rhhmm) + nminute ) * NINT(rmmss)
CALL ymds2ju( nyear, nmonth, nday, REAL(isecrst,wp), fjulday )
CALL ymds2ju( nyear, nmonth, nday, REAL(isecrst,dp), fjulday )
IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < 0.1 / rday ) fjulday = REAL(NINT(fjulday),wp) ! avoid truncation error
IF( nhour*NINT(rhhmm*rmmss) + nminute*NINT(rmmss) - ndt05 .LT. 0 ) fjulday = fjulday+1. ! move back to the day at nit000 (and not at nit000 - 1)
......@@ -115,7 +115,8 @@ CONTAINS
nday_year = nday + SUM( nmonth_len(1:nmonth - 1) )
!compute number of days between last Monday and today
CALL ymds2ju( 1900, 01, 01, 0.0_wp, zjul ) ! compute julian day value of 01.01.1900 (our reference that was a Monday)
CALL ymds2ju( 1900, 01, 01, 0.0_dp, zjul ) ! compute julian day value of 01.01.1900 (our reference that was a Monday)
inbday = FLOOR(fjulday - zjul) ! compute nb day between 01.01.1900 and start of current day
imonday = MOD(inbday, 7) ! compute nb day between last monday and current day
IF (imonday .LT. 0) imonday = imonday + 7 ! Avoid negative values for dates before 01.01.1900
......@@ -198,7 +199,7 @@ CONTAINS
nmonth_beg(jm) = nmonth_beg(jm+1) - nsecd * nmonth_len(jm)
END DO
!
END SUBROUTINE
END SUBROUTINE day_mth
SUBROUTINE day( kt )
......@@ -220,7 +221,7 @@ CONTAINS
INTEGER, INTENT(in) :: kt ! ocean time-step indices
!
CHARACTER (len=25) :: charout
REAL(wp) :: zprec ! fraction of day corresponding to 0.1 second
REAL(dp) :: zprec ! fraction of day corresponding to 0.1 second
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('day')
......@@ -259,7 +260,7 @@ CONTAINS
ndastp = nyear * 10000 + nmonth * 100 + nday ! New date
!
!compute first day of the year in julian days
CALL ymds2ju( nyear, 01, 01, 0.0_wp, fjulstartyear )
CALL ymds2ju( nyear, 01, 01, 0.0_dp, fjulstartyear )
!
IF(lwp) WRITE(numout,'(a,i8,a,i4.4,a,i2.2,a,i2.2,a,i3.3)') '======>> time-step =', kt, &
& ' New day, DATE Y/M/D = ', nyear, '/', nmonth, '/', nday, ' nday_year = ', nday_year
......@@ -310,7 +311,7 @@ CONTAINS
INTEGER , INTENT(in) :: kt ! ocean time-step
CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
!
REAL(wp) :: zkt, zndastp, zdayfrac, ksecs, ktime
REAL(dp) :: zkt, zndastp, zdayfrac, ksecs, ktime
INTEGER :: ihour, iminute, isecond
!!----------------------------------------------------------------------
......
src/OCE/DOM/depth_e3.F90
View file @ 672ff863
......@@ -143,7 +143,8 @@ CONTAINS
!!
!! ** Action : - pe3t_1d, pe3w_1d : scale factor of t- and w-point (m)
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pe3t_3d , pe3w_3d ! vert. scale factors [m]
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pe3w_3d! vert. scale factors [m]
REAL(dp), DIMENSION(:,:,:), INTENT(in ) :: pe3t_3d! vert. scale factors [m]
REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdept_3d, pdepw_3d ! depth = SUM( e3 ) [m]
!
INTEGER :: jk ! dummy loop indices
......
src/OCE/DOM/dom_oce.F90
View file @ 672ff863
......@@ -100,12 +100,16 @@ MODULE dom_oce
!! ---------------------------------------------------------------------
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: glamt , glamu, glamv , glamf !: longitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: gphit , gphiu, gphiv , gphif !: latitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t , e2t , r1_e1t, r1_e2t !: t-point horizontal scale factors [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u , e2u , r1_e1u, r1_e2u !: horizontal scale factors at u-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v , e2v , r1_e1v, r1_e2v !: horizontal scale factors at v-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f , e2f , r1_e1f, r1_e2f !: horizontal scale factors at f-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: r1_e1t, r1_e2t!: t-point horizontal scale factors [m]
REAL(dp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t, e2t!: t-point horizontal scale factors [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e2u, r1_e1u, r1_e2u!: horizontal scale factors at u-point [m]
REAL(dp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u!: horizontal scale factors at u-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v, r1_e1v, r1_e2v!: horizontal scale factors at v-point [m]
REAL(dp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e2v!: horizontal scale factors at v-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: r1_e1f, r1_e2f!: horizontal scale factors at f-point [m]
REAL(dp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f, e2f!: horizontal scale factors at f-point [m]
!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2t , r1_e1e2t !: associated metrics at t-point
REAL(dp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2t , r1_e1e2t !: associated metrics at t-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2u , r1_e1e2u , e2_e1u !: associated metrics at u-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2v , r1_e1e2v , e1_e2v !: associated metrics at v-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2f , r1_e1e2f !: associated metrics at f-point
......@@ -130,7 +134,7 @@ MODULE dom_oce
LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate
LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF
! ! reference scale factors
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0 !: t- vert. scale factor [m]
REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0 !: t- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_0 !: u- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v_0 !: v- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f_0 !: f- vert. scale factor [m]
......@@ -145,9 +149,11 @@ MODULE dom_oce
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f !: F-point vert. scale factor [m]
#endif
! ! time-dependent ratio ssh / h_0 (domqco)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: r3t, r3u, r3v !: time-dependent ratio at t-, u- and v-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: r3u, r3v!: time-dependent ratio at t-, u- and v-point [-]
REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: r3t!: time-dependent ratio at t-, u- and v-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3f !: mid-time-level ratio at f-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3t_f, r3u_f, r3v_f !: now time-filtered ratio at t-, u- and v-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3u_f, r3v_f!: now time-filtered ratio at t-, u- and v-point [-]
REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3t_f!: now time-filtered ratio at t-, u- and v-point [-]
! ! reference depths of cells
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0 !: t- depth [m]
......@@ -161,7 +167,7 @@ MODULE dom_oce
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w
#endif
! ! reference heights of ocean water column and its inverse
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0, r1_ht_0 !: t-depth [m] and [1/m]
REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0, r1_ht_0 !: t-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0, r1_hu_0 !: u-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0, r1_hv_0 !: v-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hf_0, r1_hf_0 !: f-depth [m] and [1/m]
......@@ -199,7 +205,8 @@ MODULE dom_oce
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, wmask, fmask !: land/ocean mask at T-, U-, V-, W- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wumask, wvmask !: land/ocean mask at WU- and WV-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: fe3mask !: land/ocean mask at F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_upd, umask_upd, vmask_upd !: land/ocean mask at F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_upd, umask_upd, vmask_upd !: agrif mask defining barotropic update
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_agrif !: agrif mask at T-points excluding ghosts and updated areas
!!----------------------------------------------------------------------
!! calendar variables
......@@ -215,8 +222,8 @@ MODULE dom_oce
INTEGER , PUBLIC :: nsec_month !: seconds between 00h 1st day of the current month and half of the current time step
INTEGER , PUBLIC :: nsec_monday !: seconds between 00h of the last Monday and half of the current time step
INTEGER , PUBLIC :: nsec_day !: seconds between 00h of the current day and half of the current time step
REAL(wp), PUBLIC :: fjulday !: current julian day
REAL(wp), PUBLIC :: fjulstartyear !: first day of the current year in julian days
REAL(dp), PUBLIC :: fjulday !: current julian day
REAL(dp), PUBLIC :: fjulstartyear !: first day of the current year in julian days
REAL(wp), PUBLIC :: adatrj !: number of elapsed days since the begining of the whole simulation
! !: (cumulative duration of previous runs that may have used different time-step size)
INTEGER , PUBLIC, DIMENSION( 0: 2) :: nyear_len !: length in days of the previous/current/next year
......@@ -342,7 +349,8 @@ CONTAINS
!
#if defined key_agrif
ii = ii+1
ALLOCATE( tmask_upd(jpi,jpj) , umask_upd(jpi,jpj), vmask_upd(jpi,jpj) , STAT=ierr(ii) )
ALLOCATE( tmask_upd(jpi,jpj) , umask_upd(jpi,jpj), vmask_upd(jpi,jpj), &
& tmask_agrif(jpi,jpj), STAT=ierr(ii) )
#endif
!
dom_oce_alloc = MAXVAL(ierr)
......
src/OCE/DOM/domain.F90
View file @ 672ff863
......@@ -64,6 +64,7 @@ MODULE domain
PUBLIC domain_cfg ! called by nemogcm.F90
!! * Substitutions
# include "single_precision_substitute.h90"
# include "do_loop_substitute.h90"
!!-------------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -269,10 +270,10 @@ CONTAINS
!!----------------------------------------------------------------------
!
NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list, &
& nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl , &
& nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , ln_rstdate, nn_rstctl , &
& nn_it000, nn_itend , nn_date0 , nn_time0 , nn_leapy , nn_istate , &
& nn_stock, nn_write , ln_mskland , ln_clobber , nn_chunksz, ln_1st_euler , &
& ln_cfmeta, ln_xios_read, nn_wxios
& ln_cfmeta, ln_xios_read, nn_wxios, ln_rst_eos
NAMELIST/namdom/ ln_linssh, rn_Dt, rn_atfp, ln_crs, ln_c1d, ln_meshmask
NAMELIST/namtile/ ln_tile, nn_ltile_i, nn_ltile_j
#if defined key_netcdf4
......@@ -376,9 +377,11 @@ CONTAINS
WRITE(numout,*) ' frequency of output file nn_write = ', nn_write
#endif
WRITE(numout,*) ' mask land points ln_mskland = ', ln_mskland
WRITE(numout,*) ' date-stamp restart files ln_rstdate = ', ln_rstdate
WRITE(numout,*) ' additional CF standard metadata ln_cfmeta = ', ln_cfmeta
WRITE(numout,*) ' overwrite an existing file ln_clobber = ', ln_clobber
WRITE(numout,*) ' NetCDF chunksize (bytes) nn_chunksz = ', nn_chunksz
WRITE(numout,*) ' check restart equation of state ln_rst_eos = ', ln_rst_eos
IF( TRIM(Agrif_CFixed()) == '0' ) THEN
WRITE(numout,*) ' READ restart for a single file using XIOS ln_xios_read =', ln_xios_read
WRITE(numout,*) ' Write restart using XIOS nn_wxios = ', nn_wxios
......@@ -492,6 +495,7 @@ CONTAINS
WRITE(numout,*)
IF( ln_tile ) THEN
WRITE(numout,*) ' The domain will be decomposed into tiles of size', nn_ltile_i, 'x', nn_ltile_j
CALL ctl_warn( 'dom_nam', 'you use ln_tile=.true., this may slow down NEMO performances' )
ELSE
WRITE(numout,*) ' Domain tiling will NOT be used'
ENDIF
......@@ -545,12 +549,12 @@ CONTAINS
!
llmsk = tmask_i(:,:) == 1._wp
!
CALL mpp_minloc( 'domain', glamt(:,:), llmsk, zglmin, imil )
CALL mpp_minloc( 'domain', gphit(:,:), llmsk, zgpmin, imip )
CALL mpp_minloc( 'domain', e1t(:,:), llmsk, ze1min, imi1 )
CALL mpp_minloc( 'domain', e2t(:,:), llmsk, ze2min, imi2 )
CALL mpp_maxloc( 'domain', glamt(:,:), llmsk, zglmax, imal )
CALL mpp_maxloc( 'domain', gphit(:,:), llmsk, zgpmax, imap )
CALL mpp_minloc( 'domain', CASTDP(glamt(:,:)), llmsk, zglmin, imil )
CALL mpp_minloc( 'domain', CASTDP(gphit(:,:)), llmsk, zgpmin, imip )
CALL mpp_minloc( 'domain', CASTDP(e1t(:,:)), llmsk, ze1min, imi1 )
CALL mpp_minloc( 'domain', CASTDP(e2t(:,:)), llmsk, ze2min, imi2 )
CALL mpp_maxloc( 'domain', CASTDP(glamt(:,:)), llmsk, zglmax, imal )
CALL mpp_maxloc( 'domain', CASTDP(gphit(:,:)), llmsk, zgpmax, imap )
CALL mpp_maxloc( 'domain', e1t(:,:), llmsk, ze1max, ima1 )
CALL mpp_maxloc( 'domain', e2t(:,:), llmsk, ze2max, ima2 )
!
......
src/OCE/DOM/domhgr.F90
View file @ 672ff863
......@@ -114,9 +114,9 @@ CONTAINS
! make sure that periodicities are properly applied
CALL lbc_lnk( 'dom_hgr', glamt, 'T', 1._wp, glamu, 'U', 1._wp, glamv, 'V', 1._wp, glamf, 'F', 1._wp, &
& gphit, 'T', 1._wp, gphiu, 'U', 1._wp, gphiv, 'V', 1._wp, gphif, 'F', 1._wp, &
& e1t, 'T', 1._wp, e1u, 'U', 1._wp, e1v, 'V', 1._wp, e1f, 'F', 1._wp, &
& e2t, 'T', 1._wp, e2u, 'U', 1._wp, e2v, 'V', 1._wp, e2f, 'F', 1._wp, &
& kfillmode = jpfillcopy ) ! do not put 0 over closed boundaries
& e2u, 'U', 1._wp, e1v, 'V', 1._wp, kfillmode = jpfillcopy ) ! do not put 0 over closed boundaries
CALL lbc_lnk( 'dom_hgr', e1t, 'T', 1._dp, e2t, 'T', 1._dp, e1u, 'U', 1._dp, e2v, 'V', 1._dp, &
& e1f, 'F', 1._dp, e2f, 'F', 1._dp, kfillmode = jpfillcopy)
!
ENDIF
!
......@@ -184,8 +184,10 @@ CONTAINS
REAL(wp), DIMENSION(:,:), INTENT(out) :: pphit, pphiu, pphiv, pphif ! latitude outputs
INTEGER , INTENT(out) :: kff ! =1 Coriolis parameter read here, =0 otherwise
REAL(wp), DIMENSION(:,:), INTENT(out) :: pff_f, pff_t ! Coriolis factor at f-point (if found in file)
REAL(wp), DIMENSION(:,:), INTENT(out) :: pe1t, pe1u, pe1v, pe1f ! i-scale factors
REAL(wp), DIMENSION(:,:), INTENT(out) :: pe2t, pe2u, pe2v, pe2f ! j-scale factors
REAL(wp), DIMENSION(:,:), INTENT(out) :: pe1v! i-scale factors
REAL(dp), DIMENSION(:,:), INTENT(out) :: pe1t, pe1u, pe1f! i-scale factors
REAL(wp), DIMENSION(:,:), INTENT(out) :: pe2u! j-scale factors
REAL(dp), DIMENSION(:,:), INTENT(out) :: pe2t, pe2v, pe2f! j-scale factors
INTEGER , INTENT(out) :: ke1e2u_v ! =1 u- & v-surfaces read here, =0 otherwise
REAL(wp), DIMENSION(:,:), INTENT(out) :: pe1e2u, pe1e2v ! u- & v-surfaces (if found in file)
!
......@@ -210,15 +212,15 @@ CONTAINS
CALL iom_get( inum, jpdom_global, 'gphiv', pphiv, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'gphif', pphif, cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
!
CALL iom_get( inum, jpdom_global, 'e1t' , pe1t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e1u' , pe1u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e1t' , pe1t , cd_type = 'T', psgn = 1._dp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e1u' , pe1u , cd_type = 'U', psgn = 1._dp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e1v' , pe1v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e1f' , pe1f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e1f' , pe1f , cd_type = 'F', psgn = 1._dp, kfill = jpfillcopy )
!
CALL iom_get( inum, jpdom_global, 'e2t' , pe2t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e2t' , pe2t , cd_type = 'T', psgn = 1._dp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e2u' , pe2u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e2v' , pe2v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e2f' , pe2f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e2v' , pe2v , cd_type = 'V', psgn = 1._dp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e2f' , pe2f , cd_type = 'F', psgn = 1._dp, kfill = jpfillcopy )
!
IF( iom_varid( inum, 'ff_f', ldstop = .FALSE. ) > 0 .AND. &
& iom_varid( inum, 'ff_t', ldstop = .FALSE. ) > 0 ) THEN
......
src/OCE/DOM/dommsk.F90
View file @ 672ff863
......@@ -33,6 +33,9 @@ MODULE dommsk
USE iom ! IOM library
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE lib_mpp ! Massively Parallel Processing library
USE iom ! For shlat2d
USE fldread ! for sn_shlat2d
IMPLICIT NONE
PRIVATE
......@@ -85,7 +88,11 @@ CONTAINS
INTEGER :: iktop, ikbot ! - -
INTEGER :: ios, inum
!!
NAMELIST/namlbc/ rn_shlat, ln_vorlat
REAL(wp) :: zshlat !: locally modified shlat for some strait
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zshlat2d
LOGICAL :: ln_shlat2d
CHARACTER(len = 256) :: cn_shlat2d_file, cn_shlat2d_var
NAMELIST/namlbc/ rn_shlat, ln_vorlat, ln_shlat2d, cn_shlat2d_file, cn_shlat2d_var
NAMELIST/nambdy/ ln_bdy ,nb_bdy, ln_coords_file, cn_coords_file, &
& ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, &
& cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, &
......@@ -110,12 +117,20 @@ CONTAINS
ENDIF
!
IF(lwp) WRITE(numout,*)
IF ( rn_shlat == 0. ) THEN ; IF(lwp) WRITE(numout,*) ' ==>>> ocean lateral free-slip'
ELSEIF ( rn_shlat == 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ==>>> ocean lateral no-slip'
ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ==>>> ocean lateral partial-slip'
ELSEIF ( 2. < rn_shlat ) THEN ; IF(lwp) WRITE(numout,*) ' ==>>> ocean lateral strong-slip'
ELSE
CALL ctl_stop( 'dom_msk: wrong value for rn_shlat (i.e. a negalive value). We stop.' )
IF ( ln_shlat2d ) THEN
IF(lwp) WRITE(numout,*) ' READ shlat as a 2D coefficient in a file '
ALLOCATE( zshlat2d(jpi,jpj) )
CALL iom_open(TRIM(cn_shlat2d_file), inum)
CALL iom_get (inum, jpdom_global, TRIM(cn_shlat2d_var), zshlat2d, 1) !
CALL iom_close(inum)
ELSE
IF ( rn_shlat == 0. ) THEN ; IF(lwp) WRITE(numout,*) ' ==>>> ocean lateral free-slip'
ELSEIF ( rn_shlat == 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ==>>> ocean lateral no-slip'
ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ==>>> ocean lateral partial-slip'
ELSEIF ( 2. < rn_shlat ) THEN ; IF(lwp) WRITE(numout,*) ' ==>>> ocean lateral strong-slip'
ELSE
CALL ctl_stop( 'dom_msk: wrong value for rn_shlat (i.e. a negalive value). We stop.' )
ENDIF
ENDIF
! Ocean/land mask at t-point (computed from ko_top and ko_bot)
......@@ -158,6 +173,15 @@ CONTAINS
! In case of a coarsened grid, account her for possibly aditionnal
! masked points; these have been read in the mesh file and stored in mbku, mbkv, mbkf
DO_2D( 0, 0, 0, 0 )
! Ugly patch to accomodate batropic case (jpk=2)
! but with possible masked faces in the coarsening case
! A better way would be to keep track of mbku/v/f=0 until here
! but these have been assigned a minimum of 1. TBC
IF (jpk>2) THEN
IF (mbku(ji,jj)==1) umask(ji,jj,:) = 0._wp
IF (mbkv(ji,jj)==1) vmask(ji,jj,:) = 0._wp
IF (mbkf(ji,jj)==1) fmask(ji,jj,:) = 0._wp
ENDIF
IF ( MAXVAL(umask(ji,jj,:))/=0._wp ) umask(ji,jj,mbku(ji,jj)+1:jpk) = 0._wp
IF ( MAXVAL(vmask(ji,jj,:))/=0._wp ) vmask(ji,jj,mbkv(ji,jj)+1:jpk) = 0._wp
IF ( MAXVAL(fmask(ji,jj,:))/=0._wp ) fmask(ji,jj,mbkf(ji,jj)+1:jpk) = 0._wp
......@@ -198,14 +222,26 @@ CONTAINS
! Lateral boundary conditions on velocity (modify fmask)
! ---------------------------------------
IF( rn_shlat /= 0._wp ) THEN ! Not free-slip lateral boundary condition
IF( rn_shlat /= 0._wp .or. ln_shlat2d ) THEN ! Not free-slip lateral boundary condition
!
IF ( ln_shlat2d ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
IF( fmask(ji,jj,jk) == 0._wp ) THEN
fmask(ji,jj,jk) = zshlat2d(ji,jj) * MIN( 1._wp , MAX( umask(ji,jj,jk), umask(ji,jj+1,jk), &
& vmask(ji,jj,jk), vmask(ji+1,jj,jk) ) )
ENDIF
END_3D
ELSE
DO_3D( 0, 0, 0, 0, 1, jpk )
IF( fmask(ji,jj,jk) == 0._wp ) THEN
fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( umask(ji,jj,jk), umask(ji,jj+1,jk), &
& vmask(ji,jj,jk), vmask(ji+1,jj,jk) ) )
ENDIF
END_3D
END IF
!
IF( ln_shlat2d ) DEALLOCATE( zshlat2d )
!
DO_3D( 0, 0, 0, 0, 1, jpk )
IF( fmask(ji,jj,jk) == 0._wp ) THEN
fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( umask(ji,jj,jk), umask(ji,jj+1,jk), &
& vmask(ji,jj,jk), vmask(ji+1,jj,jk) ) )
ENDIF
END_3D
CALL lbc_lnk( 'dommsk', fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
!
! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 ) depending on ln_vorlat
......@@ -215,7 +251,9 @@ CONTAINS
! User defined alteration of fmask (use to reduce ocean transport in specified straits)
! --------------------------------
!
CALL usr_def_fmask( cn_cfg, nn_cfg, fmask )
IF ( .not. ln_shlat2d ) THEN
CALL usr_def_fmask( cn_cfg, nn_cfg, fmask )
ENDIF
!
#if defined key_agrif
! Reset masks defining updated points over parent grids
......@@ -225,6 +263,10 @@ CONTAINS
tmask_upd(:,:) = 0._wp
umask_upd(:,:) = 0._wp
vmask_upd(:,:) = 0._wp
!
! Reset mask defining actual computationnal domain
! e.g. excluding ghosts and updated cells.
tmask_agrif(:,:) = 1._wp
#endif
!
END SUBROUTINE dom_msk
......
src/OCE/DOM/domqco.F90
View file @ 672ff863
......@@ -146,8 +146,9 @@ CONTAINS
!! Flux Form : simple averaging
!! - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional)
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pssh ! sea surface height [m]
REAL(wp), DIMENSION(:,:) , INTENT( out) :: pr3t, pr3u, pr3v ! ssh/h0 ratio at t-, u-, v-,points [-]
REAL(dp), DIMENSION(:,:) , INTENT(in ) :: pssh ! sea surface height [m]
REAL(wp), DIMENSION(:,:) , INTENT( out) :: pr3u, pr3v! ssh/h0 ratio at t-, u-, v-,points [-]
REAL(dp), DIMENSION(:,:) , INTENT( out) :: pr3t! ssh/h0 ratio at t-, u-, v-,points [-]
REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT( out) :: pr3f ! ssh/h0 ratio at f-point [-]
!
INTEGER :: ji, jj ! dummy loop indices
......
src/OCE/DOM/domutl.F90
View file @ 672ff863
......@@ -28,6 +28,8 @@ MODULE domutl
PUBLIC dom_uniq ! Called by dommsk and domwri
PUBLIC is_tile
!! * Substitutions
# include "single_precision_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.2 , NEMO Consortium (2020)
!! $Id: domutl.F90 14834 2021-05-11 09:24:44Z hadcv $
......@@ -75,7 +77,7 @@ CONTAINS
zgphi(:,:) = zgphi(:,:) - plat
zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:)
!
CALL mpp_minloc( 'domngb', zdist(:,:), llmsk, zmini, iloc, ldhalo = .TRUE. )
CALL mpp_minloc( 'domngb', CASTDP(zdist(:,:)), llmsk, zmini, iloc, ldhalo = .TRUE. )
kii = iloc(1)
kjj = iloc(2)
!
......@@ -107,7 +109,7 @@ CONTAINS
ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
!
puniq(:,:) = ztstref(:,:) ! default definition
CALL lbc_lnk( 'domwri', puniq, cdgrd, 1. ) ! apply boundary conditions
CALL lbc_lnk( 'domwri', puniq, cdgrd, 1._wp ) ! apply boundary conditions
lluniq(:,:,1) = puniq(:,:) == ztstref(:,:) ! check which values have not been changed
!
puniq(:,:) = REAL( COUNT( lluniq(:,:,:), dim = 3 ), wp )
......
src/OCE/DOM/domvvl.F90
View file @ 672ff863
......@@ -692,22 +692,15 @@ CONTAINS
!! - vertical interpolation: simple averaging
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3_in ! input e3 to be interpolated
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3_out ! output interpolated e3
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out) :: pe3_out ! output interpolated e3
CHARACTER(LEN=*) , INTENT(in ) :: pout ! grid point of out scale factors
! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: iku, ikum1, ikv, ikvm1, ikf, ikfm1
REAL(wp) :: zlnwd ! =1./0. when ln_wd_il = T/F
REAL(wp), DIMENSION(jpi,jpj) :: zssh ! work array to retrieve ssh (nn_vvl_interp > 1)
!!----------------------------------------------------------------------
!
IF(ln_wd_il) THEN
zlnwd = 1.0_wp
ELSE
zlnwd = 0.0_wp
END IF
!
SELECT CASE ( pout ) !== type of interpolation ==!
!
CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean
......@@ -715,7 +708,7 @@ CONTAINS
CASE ( 0 )
!
DO_3D( 1, 0, 1, 0, 1, jpk )
pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) &
pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1e2u(ji,jj) &
& * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
& + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
END_3D
......@@ -723,7 +716,7 @@ CONTAINS
CASE ( 1 )
!
DO_3D( 1, 0, 1, 0, 1, jpk )
pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) &
pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1e2u(ji,jj) &
& * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
& + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
END_3D
......@@ -732,8 +725,7 @@ CONTAINS
DO_2D( 1, 0, 1, 0 )
iku = mbku(ji ,jj)
ikum1 = iku - 1
pe3_out(ji,jj,iku) = ( umask(ji,jj,iku) * (1.0_wp - zlnwd) + zlnwd ) &
& * ( 0.5_wp * r1_e1e2u(ji,jj) &
pe3_out(ji,jj,iku) = umask(ji,jj,iku) * ( 0.5_wp * r1_e1e2u(ji,jj) &
& * ( e1e2t(ji ,jj) * ( SUM(tmask(ji ,jj,:)*(pe3_in(ji ,jj,:) - e3t_0(ji ,jj,:))) ) &
& + e1e2t(ji+1,jj) * ( SUM(tmask(ji+1,jj,:)*(pe3_in(ji+1,jj,:) - e3t_0(ji+1,jj,:))) ) ) &
& - SUM(pe3_out(ji,jj,1:ikum1)))
......@@ -743,7 +735,7 @@ CONTAINS
zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3)
!
DO_3D( 1, 0, 1, 0, 1, jpk )
pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) &
pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1e2u(ji,jj) &
& * ( e1e2t(ji ,jj) * zssh(ji ,jj) + e1e2t(ji+1,jj) * zssh(ji+1,jj)) &
& * e3u_0(ji,jj,jk) / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) )
END_3D
......@@ -758,7 +750,7 @@ CONTAINS
CASE ( 0 )
!
DO_3D( 1, 0, 1, 0, 1, jpk )
pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) &
pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1e2v(ji,jj) &
& * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
& + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
END_3D
......@@ -766,7 +758,7 @@ CONTAINS
CASE ( 1 )
!
DO_3D( 1, 0, 1, 0, 1, jpk )
pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) &
pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1e2v(ji,jj) &
& * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
& + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
END_3D
......@@ -775,8 +767,7 @@ CONTAINS
DO_2D( 1, 0, 1, 0 )
ikv = mbkv(ji ,jj)
ikvm1 = ikv - 1
pe3_out(ji,jj,ikv) = ( vmask(ji,jj,ikv) * (1.0_wp - zlnwd) + zlnwd ) &
& * ( 0.5_wp * r1_e1e2v(ji,jj) &
pe3_out(ji,jj,ikv) = vmask(ji,jj,ikv) * ( 0.5_wp * r1_e1e2v(ji,jj) &
& * ( e1e2t(ji,jj ) * ( SUM(tmask(ji,jj ,:)*(pe3_in(ji,jj ,:) - e3t_0(ji,jj ,:))) ) &
& + e1e2t(ji,jj+1) * ( SUM(tmask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3t_0(ji,jj+1,:))) ) ) &
& - SUM(pe3_out(ji,jj,1:ikvm1)))
......@@ -786,7 +777,7 @@ CONTAINS
zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3)
!
DO_3D( 1, 0, 1, 0, 1, jpk )
pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) &
pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1e2v(ji,jj) &
& * ( e1e2t(ji ,jj) * zssh(ji ,jj) + e1e2t(ji,jj+1) * zssh(ji,jj+1)) &
& * e3v_0(ji,jj,jk) / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) )
END_3D
......@@ -801,7 +792,7 @@ CONTAINS
CASE ( 0 )
!
DO_3D( 0, 0, 0, 0, 1, jpk )
pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
& * r1_e1e2f(ji,jj) &
& * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) &
& + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
......@@ -810,7 +801,7 @@ CONTAINS
CASE ( 1 )
!
DO_3D( 0, 0, 0, 0, 1, jpk )
pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
& * r1_e1e2f(ji,jj) &
& * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) &
& + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
......@@ -820,7 +811,7 @@ CONTAINS
DO_2D( 0, 0, 0, 0 )
ikf = MIN(mbku(ji ,jj),mbku(ji,jj+1))
ikfm1 = ikf - 1
pe3_out(ji,jj,ikf) = ( umask(ji,jj,ikf) * umask(ji,jj+1,ikf) * (1.0_wp - zlnwd) + zlnwd ) &
pe3_out(ji,jj,ikf) = umask(ji,jj,ikf) * umask(ji,jj+1,ikf) &
& * ( 0.5_wp * r1_e1e2f(ji,jj) &
& * ( e1e2u(ji,jj ) * ( SUM(umask(ji,jj ,:)*(pe3_in(ji,jj ,:) - e3u_0(ji,jj ,:))) ) &
& + e1e2u(ji,jj+1) * ( SUM(umask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3u_0(ji,jj+1,:))) ) ) &
......@@ -831,7 +822,7 @@ CONTAINS
zssh(:,:) = SUM(umask(:,:,:)*(pe3_in(:,:,:)-e3u_0(:,:,:)), DIM=3)
!
DO_3D( 0, 0, 0, 0, 1, jpk )
pe3_out(ji,jj,jk) = ( umask(ji,jj,jk)* umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
pe3_out(ji,jj,jk) = umask(ji,jj,jk)* umask(ji,jj+1,jk) &
& * 0.5_wp * r1_e1e2f(ji,jj) &
& * (e1e2u(ji ,jj) * zssh(ji ,jj) + e1e2u(ji,jj+1) * zssh(ji,jj+1)) &
& * e3f_0(ji,jj,jk) / ( hf_0(ji,jj) + 1._wp - ssumask(ji,jj)*ssumask(ji,jj+1) )
......@@ -848,9 +839,9 @@ CONTAINS
! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
!!gm BUG? use here wmask in case of ISF ? to be checked
DO jk = 2, jpk
pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
& * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) &
& + 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) &
& * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) &
& + 0.5_wp * tmask(:,:,jk) &
& * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) )
END DO
!
......@@ -860,9 +851,9 @@ CONTAINS
! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
!!gm BUG? use here wumask in case of ISF ? to be checked
DO jk = 2, jpk
pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
& * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) &
& + 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) &
& * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) &
& + 0.5_wp * umask(:,:,jk) &
& * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) )
END DO
!
......@@ -872,9 +863,9 @@ CONTAINS
! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
!!gm BUG? use here wvmask in case of ISF ? to be checked
DO jk = 2, jpk
pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
& * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) &
& + 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) &
& * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) &
& + 0.5_wp * vmask(:,:,jk) &
& * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) )
END DO
END SELECT
......
src/OCE/DOM/domzgr.F90
View file @ 672ff863
......@@ -113,10 +113,11 @@ CONTAINS
& k_top , k_bot ) ! 1st & last ocean level
!
! make sure that periodicities are properly applied
CALL lbc_lnk( 'dom_zgr', gdept_0, 'T', 1._wp, gdepw_0, 'W', 1._wp, &
& e3t_0, 'T', 1._wp, e3u_0, 'U', 1._wp, e3v_0, 'V', 1._wp, e3f_0, 'F', 1._wp, &
CALL lbc_lnk( 'dom_zgr', gdept_0, 'T', 1._wp, gdepw_0, 'W', 1._wp, &
& e3u_0, 'U', 1._wp, e3v_0, 'V', 1._wp, e3f_0, 'F', 1._wp, &
& e3w_0, 'W', 1._wp, e3uw_0, 'U', 1._wp, e3vw_0, 'V', 1._wp, &
& kfillmode = jpfillcopy ) ! do not put 0 over closed boundaries
CALL lbc_lnk( 'dom_zgr', e3t_0, 'T', 1._dp, kfillmode = jpfillcopy )
ztopbot(:,:,1) = REAL(k_top, wp)
ztopbot(:,:,2) = REAL(k_bot, wp)
CALL lbc_lnk( 'dom_zgr', ztopbot, 'T', 1._wp, kfillmode = jpfillcopy ) ! do not put 0 over closed boundaries
......@@ -139,7 +140,7 @@ CONTAINS
IF( .NOT. ( l_Jperio .OR. l_NFold ) ) THEN ! N closed:
zmsk(:,mj0(jpjglo-nn_hls):mj1(jpjglo-nn_hls) ) = 0._wp ! last line of inner global domain at 0
ENDIF
CALL lbc_lnk( 'usrdef_zgr', zmsk, 'T', 1. ) ! set halos
CALL lbc_lnk( 'usrdef_zgr', zmsk, 'T', 1._wp ) ! set halos
k_top(:,:) = k_top(:,:) * NINT( zmsk(:,:) )
!
!!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears
......@@ -235,7 +236,8 @@ CONTAINS
REAL(wp), DIMENSION(:) , INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m]
REAL(wp), DIMENSION(:) , INTENT(out) :: pe3t_1d , pe3w_1d ! 1D vertical scale factors [m]
REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pdept, pdepw ! grid-point depth [m]
REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3t , pe3u , pe3v , pe3f ! vertical scale factors [m]
REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3u, pe3v, pe3f! vertical scale factors [m]
REAL(dp), DIMENSION(:,:,:), INTENT(out) :: pe3t! vertical scale factors [m]
REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3w , pe3uw, pe3vw ! - - -
INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top , k_bot ! first & last ocean level
INTEGER , INTENT(out) :: k_mbkuvf ! ==1 if mbku, mbkv, mbkf are in file
......@@ -286,7 +288,7 @@ CONTAINS
CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , pe3t_1d ) ! 1D reference coordinate
CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , pe3w_1d )
!
CALL iom_get( inum, jpdom_global, 'e3t_0' , pe3t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy ) ! 3D coordinate
CALL iom_get( inum, jpdom_global, 'e3t_0' , pe3t , cd_type = 'T', psgn = 1._dp, kfill = jpfillcopy ) ! 3D coordinate
CALL iom_get( inum, jpdom_global, 'e3u_0' , pe3u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e3v_0' , pe3v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e3f_0' , pe3f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
......
src/OCE/DOM/domzgr_substitute.h90
View file @ 672ff863
......@@ -51,4 +51,3 @@
# define gde3w(i,j,k) gdept_0(i,j,k)
#endif
!!----------------------------------------------------------------------
src/OCE/DOM/dtatsd.F90
View file @ 672ff863
......@@ -135,13 +135,13 @@ CONTAINS
!! ** Action : ptsd T-S data on medl mesh and interpolated at time-step kt
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step
REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts), INTENT( out) :: ptsd ! T & S data
REAL(dp), DIMENSION(A2D(nn_hls),jpk,jpts), INTENT( out) :: ptsd ! T & S data
!
INTEGER :: ji, jj, jk, jl, jkk ! dummy loop indicies
INTEGER :: ik, il0, il1, ii0, ii1, ij0, ij1 ! local integers
INTEGER, DIMENSION(jpts), SAVE :: irec_b, irec_n
REAL(wp):: zl, zi ! local scalars
REAL(wp), DIMENSION(jpk) :: ztp, zsp ! 1D workspace
REAL(dp):: zl, zi ! local scalars
REAL(dp), DIMENSION(jpk) :: ztp, zsp ! 1D workspace
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only for the full domain
......
src/OCE/DOM/istate.F90
View file @ 672ff863
......@@ -132,8 +132,8 @@ CONTAINS
END DO
CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb) )
! make sure that periodicities are properly applied
CALL lbc_lnk( 'istate', ts(:,:,:,jp_tem,Kbb), 'T', 1._wp, ts(:,:,:,jp_sal,Kbb), 'T', 1._wp, &
& uu(:,:,:, Kbb), 'U', -1._wp, vv(:,:,:, Kbb), 'V', -1._wp )
CALL lbc_lnk( 'istate', ts(:,:,:,jp_tem,Kbb), 'T', 1._dp, ts(:,:,:,jp_sal,Kbb), 'T', 1._dp, &
& uu(:,:,:, Kbb), 'U', -1._dp, vv(:,:,:, Kbb), 'V', -1._dp )
ENDIF
ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones
uu (:,:,:,Kmm) = uu (:,:,:,Kbb)
......
src/OCE/DYN/dynadv.F90
View file @ 672ff863
......@@ -16,6 +16,7 @@ MODULE dynadv
USE dom_oce ! ocean space and time domain
USE dynadv_cen2 ! centred flux form advection (dyn_adv_cen2 routine)
USE dynadv_ubs ! UBS flux form advection (dyn_adv_ubs routine)
USE dynadv_up3 ! UBS flux form advection (NEW) (dyn_adv_up3 routine)
USE dynkeg ! kinetic energy gradient (dyn_keg routine)
USE dynzad ! vertical advection (dyn_zad routine)
!
......@@ -34,14 +35,16 @@ MODULE dynadv
LOGICAL, PUBLIC :: ln_dynadv_vec !: vector form
INTEGER, PUBLIC :: nn_dynkeg !: scheme of grad(KE): =0 C2 ; =1 Hollingsworth
LOGICAL, PUBLIC :: ln_dynadv_cen2 !: flux form - 2nd order centered scheme flag
LOGICAL, PUBLIC :: ln_dynadv_ubs !: flux form - 3rd order UBS scheme flag
LOGICAL, PUBLIC :: ln_dynadv_ubs !: flux form - 3rd order UBS scheme flag (OLD)
LOGICAL, PUBLIC :: ln_dynadv_up3 !: flux form - 3rd order UBS scheme flag (NEW)
INTEGER, PUBLIC :: n_dynadv !: choice of the formulation and scheme for momentum advection
! ! associated indices:
INTEGER, PUBLIC, PARAMETER :: np_LIN_dyn = 0 ! no advection: linear dynamics
INTEGER, PUBLIC, PARAMETER :: np_VEC_c2 = 1 ! vector form : 2nd order centered scheme
INTEGER, PUBLIC, PARAMETER :: np_FLX_c2 = 2 ! flux form : 2nd order centered scheme
INTEGER, PUBLIC, PARAMETER :: np_FLX_ubs = 3 ! flux form : 3rd order Upstream Biased Scheme
INTEGER, PUBLIC, PARAMETER :: np_FLX_ubs = 3 ! flux form : 3rd order Upstream Biased Scheme (OLD)
INTEGER, PUBLIC, PARAMETER :: np_FLX_up3 = 4 ! flux form : 3rd order Upstream Biased Scheme (NEW)
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -65,7 +68,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start( 'dyn_adv' )
......@@ -77,7 +80,9 @@ CONTAINS
CASE( np_FLX_c2 )
CALL dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs ) ! 2nd order centered scheme
CASE( np_FLX_ubs )
CALL dyn_adv_ubs ( kt, Kbb, Kmm, puu, pvv, Krhs ) ! 3rd order UBS scheme (UP3)
CALL dyn_adv_ubs ( kt, Kbb, Kmm, puu, pvv, Krhs ) ! 3rd order OLD UBS scheme (UP3)
CASE( np_FLX_up3 )
CALL dyn_adv_up3 ( kt, Kbb, Kmm, puu, pvv, Krhs ) ! 3rd order NEW UBS scheme (UP3)
END SELECT
!
IF( ln_timing ) CALL timing_stop( 'dyn_adv' )
......@@ -94,7 +99,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER :: ioptio, ios ! Local integer
!
NAMELIST/namdyn_adv/ ln_dynadv_OFF, ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2, ln_dynadv_ubs
NAMELIST/namdyn_adv/ ln_dynadv_OFF, ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2, ln_dynadv_ubs, ln_dynadv_up3
!!----------------------------------------------------------------------
!
IF(lwp) THEN
......@@ -115,7 +120,8 @@ CONTAINS
WRITE(numout,*) ' Vector form: 2nd order centered scheme ln_dynadv_vec = ', ln_dynadv_vec
WRITE(numout,*) ' with Hollingsworth scheme (=1) or not (=0) nn_dynkeg = ', nn_dynkeg
WRITE(numout,*) ' flux form: 2nd order centred scheme ln_dynadv_cen2 = ', ln_dynadv_cen2
WRITE(numout,*) ' 3rd order UBS scheme ln_dynadv_ubs = ', ln_dynadv_ubs
WRITE(numout,*) ' 3rd order UBS scheme (OLD) ln_dynadv_ubs = ', ln_dynadv_ubs
WRITE(numout,*) ' 3rd order UBS scheme (NEW) ln_dynadv_up3 = ', ln_dynadv_up3
ENDIF
ioptio = 0 ! parameter control and set n_dynadv
......@@ -123,6 +129,7 @@ CONTAINS
IF( ln_dynadv_vec ) THEN ; ioptio = ioptio + 1 ; n_dynadv = np_VEC_c2 ; ENDIF
IF( ln_dynadv_cen2 ) THEN ; ioptio = ioptio + 1 ; n_dynadv = np_FLX_c2 ; ENDIF
IF( ln_dynadv_ubs ) THEN ; ioptio = ioptio + 1 ; n_dynadv = np_FLX_ubs ; ENDIF
IF( ln_dynadv_up3 ) THEN ; ioptio = ioptio + 1 ; n_dynadv = np_FLX_up3 ; ENDIF
IF( ioptio /= 1 ) CALL ctl_stop( 'choose ONE and only ONE advection scheme' )
IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW ) CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' )
......@@ -138,7 +145,8 @@ CONTAINS
IF( nn_dynkeg == nkeg_C2 ) WRITE(numout,*) ' with Centered standard keg scheme'
IF( nn_dynkeg == nkeg_HW ) WRITE(numout,*) ' with Hollingsworth keg scheme'
CASE( np_FLX_c2 ) ; WRITE(numout,*) ' ==>>> flux form : 2nd order scheme is used'
CASE( np_FLX_ubs ) ; WRITE(numout,*) ' ==>>> flux form : UBS scheme is used'
CASE( np_FLX_ubs ) ; WRITE(numout,*) ' ==>>> flux form : OLD UBS scheme is used'
CASE( np_FLX_up3 ) ; WRITE(numout,*) ' ==>>> flux form : NEW UBS scheme is used'
END SELECT
ENDIF
!
......
src/OCE/DYN/dynadv_cen2.F90
View file @ 672ff863
......@@ -48,11 +48,13 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t, zfu_f, zfu_uw, zfu
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t, zfv_f, zfv_vw, zfv, zfw
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_f, zfu
REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t, zfu_uw
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_f, zfv, zfw
REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t, zfv_vw
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......
src/OCE/DYN/dynadv_ubs.F90
View file @ 672ff863
......@@ -70,12 +70,14 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t, zfu_f, zfu_uw, zfu
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t, zfv_f, zfv_vw, zfv, zfw
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_f, zfu
REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t, zfu_uw
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_f, zfv, zfw
REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t, zfv_vw
REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: zlu_uu, zlu_uv
REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: zlv_vv, zlv_vu
!!----------------------------------------------------------------------
......
src/OCE/DYN/dynadv_up3.F90 0 → 100644
View file @ 672ff863
MODULE dynadv_up3
!!======================================================================
!! *** MODULE dynadv_up3 ***
!! Ocean dynamics: Update the momentum trend with the flux form advection
!! trend using a 3rd order upstream biased scheme
!!======================================================================
!! History : 2.0 ! 2006-08 (R. Benshila, L. Debreu) Original code
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! dyn_adv_up3 : flux form momentum advection using (ln_dynadv=T)
!! an 3rd order Upstream Biased Scheme or Quick scheme
!! combined with 2nd or 4th order finite differences
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
USE trd_oce ! trends: ocean variables
USE trddyn ! trend manager: dynamics
!
USE in_out_manager ! I/O manager
USE prtctl ! Print control
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE lib_mpp ! MPP library
IMPLICIT NONE
PRIVATE
REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp ! =1/4 quick ; =1/3 3rd order UBS
PUBLIC dyn_adv_up3 ! routine called by step.F90
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynadv_up3.F90 14834 2021-05-11 09:24:44Z hadcv $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv_up3( kt, Kbb, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_up3 ***
!!
!! ** Purpose : Compute the now momentum advection trend in flux form
!! and the general trend of the momentum equation.
!!
!! ** Method : The scheme is the one implemeted in ROMS. It depends
!! on two parameter gamma1 and gamma2. The former control the
!! upstream baised part of the scheme and the later the centred
!! part: gamma1 = 0 pure centered (no diffusive part)
!! = 1/4 Quick scheme
!! = 1/3 3rd order Upstream biased scheme
!! For stability reasons, the first term of the fluxes which cor-
!! responds to a second order centered scheme is evaluated using
!! the now velocity (centered in time) while the second term which
!! is the diffusive part of the scheme, is evaluated using the
!! before velocity (forward in time).
!! Default value (hard coded in the begining of the module) are
!! gamma1=1/3 and gamma2=1/32.
!!
!! ** Action : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the 3D advective momentum trends
!!
!! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices
REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! local scalars
REAL(wp) :: zfwi, zzfu_uw_kp1, zlu_uw_kp1 ! local scalars
REAL(wp) :: zfwj, zzfv_vw_kp1, zlv_vw_kp1 ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_f, zfu
REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_f, zfv
REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t
REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: ztrdu, ztrdv
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zlu_uu, zlu_uv
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zlv_vv, zlv_vu
REAL(dp), DIMENSION(A2D(nn_hls)) :: zfw, zfu_uw, zfv_vw, zlu_uw, zlv_vw
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
!
zfu_t(:,:,:) = 0._wp
zfv_t(:,:,:) = 0._wp
zfu_f(:,:,:) = 0._wp
zfv_f(:,:,:) = 0._wp
!
zlu_uu(:,:,:) = 0._wp
zlv_vv(:,:,:) = 0._wp
zlu_uv(:,:,:) = 0._wp
zlv_vu(:,:,:) = 0._wp
!
IF( l_trddyn ) THEN ! trends: store the input trends
ztrdu(:,:,:) = puu(:,:,:,Krhs)
ztrdv(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
! ! =========================== !
DO jk = 1, jpkm1 ! Laplacian of the velocity !
! ! =========================== !
! ! horizontal volume fluxes
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zfu(ji,jj,jk) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
zfv(ji,jj,jk) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
END_2D
!
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! laplacian
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zlu_uu(ji,jj,jk) = ( ( puu (ji+1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& + ( puu (ji-1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& ) * umask(ji ,jj ,jk)
zlv_vv(ji,jj,jk) = ( ( pvv (ji ,jj+1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& + ( pvv (ji ,jj-1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& ) * vmask(ji ,jj ,jk)
zlu_uv(ji,jj,jk) = ( puu (ji ,jj+1,jk,Kbb) - puu (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
& - ( puu (ji ,jj ,jk,Kbb) - puu (ji ,jj-1,jk,Kbb) ) * fmask(ji ,jj-1,jk)
zlv_vu(ji,jj,jk) = ( pvv (ji+1,jj ,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
& - ( pvv (ji ,jj ,jk,Kbb) - pvv (ji-1,jj ,jk,Kbb) ) * fmask(ji-1,jj ,jk)
END_2D
END DO
!
IF( nn_hls == 1 ) CALL lbc_lnk( 'dynadv_ubs', zlu_uu(:,:,:), 'U', -1.0_wp , zlu_uv(:,:,:), 'U', -1.0_wp, &
& zlv_vv(:,:,:), 'V', -1.0_wp , zlv_vu(:,:,:), 'V', -1.0_wp )
!
! ! ====================== !
! ! Horizontal advection !
DO jk = 1, jpkm1 ! ====================== !
! ! horizontal volume fluxes
DO_2D( 1, 1, 1, 1 )
zfu(ji,jj,jk) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
END_2D
!
DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes at T- and F-point
zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
!
IF( zui > 0 ) THEN ; zl_u = zlu_uu(ji ,jj,jk)
ELSE ; zl_u = zlu_uu(ji+1,jj,jk)
ENDIF
IF( zvj > 0 ) THEN ; zl_v = zlv_vv(ji,jj ,jk)
ELSE ; zl_v = zlv_vv(ji,jj+1,jk)
ENDIF
!
zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj ,jk) ) &
& * ( zui - gamma1 * zl_u)
zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji ,jj+1,jk) ) &
& * ( zvj - gamma1 * zl_v)
!
zfuj = ( zfu(ji,jj,jk) + zfu(ji ,jj+1,jk) )
zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj ,jk) )
IF( zfuj > 0 ) THEN ; zl_v = zlv_vu( ji ,jj ,jk )
ELSE ; zl_v = zlv_vu( ji+1,jj ,jk )
ENDIF
IF( zfvi > 0 ) THEN ; zl_u = zlu_uv( ji,jj ,jk )
ELSE ; zl_u = zlu_uv( ji,jj+1,jk )
ENDIF
!
zfv_f(ji ,jj ,jk) = ( zfvi ) &
& * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) - gamma1 * zl_u )
zfu_f(ji ,jj ,jk) = ( zfuj ) &
& * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) - gamma1 * zl_v )
END_2D
DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) &
& + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) &
& + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
END DO
IF( l_trddyn ) THEN ! trends: send trends to trddyn for diagnostic
ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt, Kmm )
zfu_t(:,:,:) = puu(:,:,:,Krhs)
zfv_t(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
! ! ==================== !
! ! Vertical advection !
! ! ==================== !
!
! != surface vertical fluxes =! (jk = 1)
!
DO_2D( 0, 1, 0, 1 )
zfw(ji,jj) = e1e2t(ji,jj) * ww(ji,jj,1)
END_2D
!
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0.5_wp * ( zfw(ji,jj) + zfw(ji+1,jj ) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj) = 0.5_wp * ( zfw(ji,jj) + zfw(ji ,jj+1) ) * pvv(ji,jj,1,Kmm)
END_2D
ELSE ! non linear free: surface advective fluxes set to zero
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0._wp
zfv_vw(ji,jj) = 0._wp
END_2D
ENDIF
!
DO_2D( 0, 0, 0, 0 ) ! uniform shear in the 1st layer : dz(u(k=1)) = dz(u(k=2)) ==> zlu = 0
zlu_uw(ji,jj) = 0._wp
zlv_vw(ji,jj) = 0._wp
END_2D
!
DO jk = 1, jpk-2 != divergence of advective fluxes =! (jk = 1 to jpk-2)
DO_2D( 0, 1, 0, 1 ) ! vertical transport at level k+1
zfw(ji,jj) = e1e2t(ji,jj) * ww(ji,jj,jk+1)
END_2D
!
DO_2D( 0, 0, 0, 0 )
zlu_uw_kp1 = ( puu(ji,jj,jk ,Kbb) - puu(ji,jj,jk+1,Kbb) ) * wumask(ji,jj,jk+1) &
& - ( puu(ji,jj,jk+1,Kbb) - puu(ji,jj,jk+2,Kbb) ) * wumask(ji,jj,jk+2)
zlv_vw_kp1 = ( pvv(ji,jj,jk ,Kbb) - pvv(ji,jj,jk+1,Kbb) ) * wvmask(ji,jj,jk+1) &
& - ( pvv(ji,jj,jk+1,Kbb) - pvv(ji,jj,jk+2,Kbb) ) * wvmask(ji,jj,jk+2)
!
zfwi = zfw(ji,jj) + zfw(ji+1,jj)
IF( zfwi > 0 ) THEN ; zl_u = zlu_uw_kp1
ELSE ; zl_u = zlu_uw(ji,jj)
ENDIF
zfwj = zfw(ji,jj) + zfw(ji,jj+1)
IF( zfwj > 0 ) THEN ; zl_v = zlv_vw_kp1
ELSE ; zl_v = zlv_vw(ji,jj)
ENDIF
! ! vertical flux at level k+1
zzfu_uw_kp1 = 0.25_wp * ( zfw(ji,jj) + zfw(ji+1,jj ) ) * ( puu(ji,jj,jk+1,Kmm) + puu(ji,jj,jk,Kmm) - gamma1 * zl_u )
zzfv_vw_kp1 = 0.25_wp * ( zfw(ji,jj) + zfw(ji ,jj+1) ) * ( pvv(ji,jj,jk+1,Kmm) + pvv(ji,jj,jk,Kmm) - gamma1 * zl_v )
! ! divergence of vertical momentum flux
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj) - zzfu_uw_kp1 ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj) - zzfv_vw_kp1 ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
! ! store vertical flux for next level calculation
zfu_uw(ji,jj) = zzfu_uw_kp1
zfv_vw(ji,jj) = zzfv_vw_kp1
!
zlu_uw(ji,jj) = zlu_uw_kp1
zlv_vw(ji,jj) = zlv_vw_kp1
END_2D
END DO
!
jk = jpkm1 != compute last level =! (zzfu_uw_kp1 = zzfu_uw_kp1 = 0)
DO_2D( 0, 0, 0, 0 )
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - zfu_uw(ji,jj) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - zfv_vw(ji,jj) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
END_2D
!
IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic
zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
END SUBROUTINE dyn_adv_up3
!!==============================================================================
END MODULE dynadv_up3
src/OCE/DYN/dynatf.F90
View file @ 672ff863
......@@ -80,6 +80,7 @@ CONTAINS
!! * Substitutions
# include "do_loop_substitute.h90"
# include "single_precision_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynatf.F90 14834 2021-05-11 09:24:44Z hadcv $
......@@ -354,8 +355,8 @@ CONTAINS
ENDIF
ENDIF
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt - puu(:,:,:,Kaa): ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): ' , mask2=vmask )
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=CASTDP(puu(:,:,:,Kaa)), clinfo1=' nxt - puu(:,:,:,Kaa): ', mask1=umask, &
& tab3d_2=CASTDP(pvv(:,:,:,Kaa)), clinfo2=' pvv(:,:,:,Kaa): ' , mask2=vmask )
!
IF( ln_dynspg_ts ) DEALLOCATE( zue, zve )
IF( l_trddyn ) DEALLOCATE( zua, zva )
......
src/OCE/DYN/dynatf_qco.F90
View file @ 672ff863
......@@ -94,13 +94,13 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! before and after time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities to be time filtered
REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities to be time filtered
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zue3a, zue3n, zue3b, zcoef ! local scalars
REAL(wp) :: zve3a, zve3n, zve3b ! - -
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zua, zva
REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: zua, zva
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zutau, zvtau
!!----------------------------------------------------------------------
!
......@@ -195,7 +195,7 @@ CONTAINS
ENDIF ! .NOT. l_1st_euler
!
! This is needed for dyn_ldf_blp to be restartable
IF( nn_hls == 2 ) CALL lbc_lnk( 'dynatfqco', puu(:,:,:,Kmm), 'U', -1.0_wp, pvv(:,:,:,Kmm), 'V', -1.0_wp )
IF( nn_hls == 2 ) CALL lbc_lnk( 'dynatfqco', puu(:,:,:,Kmm), 'U', -1.0_dp, pvv(:,:,:,Kmm), 'V', -1.0_dp )
! Set "now" and "before" barotropic velocities for next time step:
! JC: Would be more clever to swap variables than to make a full vertical
......