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 1401 additions and 131 deletions
src/ICE/iceupdate.F90
View file @ 4a9f4b18
......@@ -383,8 +383,8 @@ CONTAINS
zat_v = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji ,jj+1 ) * tmask(ji ,jj+1,1) ) &
& / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji ,jj+1,1) )
! ! linearized quadratic drag formulation
zutau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
zvtau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
zutau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - zflagi * pu_oce(ji,jj) )
zvtau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - zflagi * pv_oce(ji,jj) )
! ! stresses at the ocean surface
utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
......
src/ICE/icewri.F90
View file @ 4a9f4b18
......@@ -103,7 +103,7 @@ CONTAINS
IF( iom_use('icethic' ) ) CALL iom_put( 'icethic', hm_i * zmsk00 ) ! ice thickness
IF( iom_use('snwthic' ) ) CALL iom_put( 'snwthic', hm_s * zmsk00 ) ! snw thickness
IF( iom_use('icebrv' ) ) CALL iom_put( 'icebrv' , bvm_i* 100. * zmsk00 ) ! brine volume
IF( iom_use('iceage' ) ) CALL iom_put( 'iceage' , om_i / rday * zmsk15 + zmiss_val * ( 1._wp - zmsk15 ) ) ! ice age
IF( iom_use('iceage' ) ) CALL iom_put( 'iceage' , om_i / rday * zmsk15 ) ! ice age
IF( iom_use('icehnew' ) ) CALL iom_put( 'icehnew', ht_i_new ) ! new ice thickness formed in the leads
IF( iom_use('snwvolu' ) ) CALL iom_put( 'snwvolu', vt_s * zmsksn ) ! snow volume
IF( iom_use('icefrb' ) ) THEN ! Ice freeboard
......@@ -118,14 +118,15 @@ CONTAINS
IF( iom_use('icehlid' ) ) CALL iom_put( 'icehlid', hm_il * zmsk00 ) ! melt pond lid depth
IF( iom_use('icevlid' ) ) CALL iom_put( 'icevlid', vt_il * zmsk00 ) ! melt pond lid total volume per unit area
! salt
IF( iom_use('icesalt' ) ) CALL iom_put( 'icesalt', sm_i * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! mean ice salinity
IF( iom_use('icesalt' ) ) CALL iom_put( 'icesalt', sm_i * zmsk00 ) ! mean ice salinity
IF( iom_use('icesalm' ) ) CALL iom_put( 'icesalm', st_i * rhoi * 1.0e-3 * zmsk00 ) ! Mass of salt in sea ice per cell area
IF( iom_use('iceepnd' ) ) CALL iom_put( 'iceepnd', SUM( a_ip_eff * a_i, dim=3 ) * zmsk00 ) ! melt pond total effective fraction per cell area
! heat
IF( iom_use('icetemp' ) ) CALL iom_put( 'icetemp', ( tm_i - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! ice mean temperature
IF( iom_use('snwtemp' ) ) CALL iom_put( 'snwtemp', ( tm_s - rt0 ) * zmsksn + zmiss_val * ( 1._wp - zmsksn ) ) ! snw mean temperature
IF( iom_use('icettop' ) ) CALL iom_put( 'icettop', ( tm_su - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! temperature at the ice surface
IF( iom_use('icetbot' ) ) CALL iom_put( 'icetbot', ( t_bo - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! temperature at the ice bottom
IF( iom_use('icetsni' ) ) CALL iom_put( 'icetsni', ( tm_si - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! temperature at the snow-ice interface
IF( iom_use('icetemp' ) ) CALL iom_put( 'icetemp', ( tm_i - rt0 ) * zmsk00 ) ! ice mean temperature
IF( iom_use('snwtemp' ) ) CALL iom_put( 'snwtemp', ( tm_s - rt0 ) * zmsksn ) ! snw mean temperature
IF( iom_use('icettop' ) ) CALL iom_put( 'icettop', ( tm_su - rt0 ) * zmsk00 ) ! temperature at the ice surface
IF( iom_use('icetbot' ) ) CALL iom_put( 'icetbot', ( t_bo - rt0 ) * zmsk00 ) ! temperature at the ice bottom
IF( iom_use('icetsni' ) ) CALL iom_put( 'icetsni', ( tm_si - rt0 ) * zmsk00 ) ! temperature at the snow-ice interface
IF( iom_use('icehc' ) ) CALL iom_put( 'icehc' , -et_i * zmsk00 ) ! ice heat content
IF( iom_use('snwhc' ) ) CALL iom_put( 'snwhc' , -et_s * zmsksn ) ! snow heat content
! momentum
......@@ -169,16 +170,16 @@ CONTAINS
! --- category-dependent fields --- !
IF( iom_use('icemask_cat' ) ) CALL iom_put( 'icemask_cat' , zmsk00l ) ! ice mask 0%
IF( iom_use('iceconc_cat' ) ) CALL iom_put( 'iceconc_cat' , a_i * zmsk00l ) ! area for categories
IF( iom_use('icethic_cat' ) ) CALL iom_put( 'icethic_cat' , h_i * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! thickness for categories
IF( iom_use('snwthic_cat' ) ) CALL iom_put( 'snwthic_cat' , h_s * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow depth for categories
IF( iom_use('icesalt_cat' ) ) CALL iom_put( 'icesalt_cat' , s_i * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! salinity for categories
IF( iom_use('iceage_cat' ) ) CALL iom_put( 'iceage_cat' , o_i / rday * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice age
IF( iom_use('icethic_cat' ) ) CALL iom_put( 'icethic_cat' , h_i * zmsk00l ) ! thickness for categories
IF( iom_use('snwthic_cat' ) ) CALL iom_put( 'snwthic_cat' , h_s * zmsksnl ) ! snow depth for categories
IF( iom_use('icesalt_cat' ) ) CALL iom_put( 'icesalt_cat' , s_i * zmsk00l ) ! salinity for categories
IF( iom_use('iceage_cat' ) ) CALL iom_put( 'iceage_cat' , o_i / rday * zmsk00l ) ! ice age
IF( iom_use('icetemp_cat' ) ) CALL iom_put( 'icetemp_cat' , ( SUM( t_i, dim=3 ) * r1_nlay_i - rt0 ) &
& * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice temperature
& * zmsk00l ) ! ice temperature
IF( iom_use('snwtemp_cat' ) ) CALL iom_put( 'snwtemp_cat' , ( SUM( t_s, dim=3 ) * r1_nlay_s - rt0 ) &
& * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow temperature
IF( iom_use('icettop_cat' ) ) CALL iom_put( 'icettop_cat' , ( t_su - rt0 ) * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! surface temperature
IF( iom_use('icebrv_cat' ) ) CALL iom_put( 'icebrv_cat' , bv_i * 100. * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume
& * zmsksnl ) ! snow temperature
IF( iom_use('icettop_cat' ) ) CALL iom_put( 'icettop_cat' , ( t_su - rt0 ) * zmsk00l ) ! surface temperature
IF( iom_use('icebrv_cat' ) ) CALL iom_put( 'icebrv_cat' , bv_i * 100. * zmsk00l ) ! brine volume
IF( iom_use('iceapnd_cat' ) ) CALL iom_put( 'iceapnd_cat' , a_ip * zmsk00l ) ! melt pond frac for categories
IF( iom_use('icevpnd_cat' ) ) CALL iom_put( 'icevpnd_cat' , v_ip * zmsk00l ) ! melt pond volume for categories
IF( iom_use('icehpnd_cat' ) ) CALL iom_put( 'icehpnd_cat' , h_ip * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond thickness for categories
......
src/OCE/ASM/asminc.F90
View file @ 4a9f4b18
......@@ -372,8 +372,8 @@ CONTAINS
ENDIF
IF ( ln_dyninc ) THEN
CALL iom_get( inum, jpdom_auto, 'bckinu', u_bkginc, 1 )
CALL iom_get( inum, jpdom_auto, 'bckinv', v_bkginc, 1 )
CALL iom_get( inum, jpdom_auto, 'bckinu', u_bkginc, 1, cd_type = 'U', psgn = -1._wp )
CALL iom_get( inum, jpdom_auto, 'bckinv', v_bkginc, 1, cd_type = 'V', psgn = -1._wp )
! Apply the masks
u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
......@@ -474,8 +474,8 @@ CONTAINS
ENDIF
!
IF ( ln_dyninc ) THEN
CALL iom_get( inum, jpdom_auto, 'un', u_bkg, cd_type = 'U', psgn = 1._wp )
CALL iom_get( inum, jpdom_auto, 'vn', v_bkg, cd_type = 'V', psgn = 1._wp )
CALL iom_get( inum, jpdom_auto, 'un', u_bkg, cd_type = 'U', psgn = -1._wp )
CALL iom_get( inum, jpdom_auto, 'vn', v_bkg, cd_type = 'V', psgn = -1._wp )
u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
ENDIF
......
src/OCE/DIA/diahth.F90
View file @ 4a9f4b18
......@@ -15,6 +15,7 @@ MODULE diahth
USE oce ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
USE phycst ! physical constants
USE zdfmxl, ONLY: zdf_mxl_zint
!
USE in_out_manager ! I/O manager
USE lib_mpp ! MPP library
......@@ -292,6 +293,9 @@ CONTAINS
!
ENDIF
! Vertically-interpolated mixed-layer depth diagnostic
CALL zdf_mxl_zint( kt, Kmm )
!
IF( ln_timing ) CALL timing_stop('dia_hth')
!
......
src/OCE/DIA/diaprod.F90 0 → 100644
View file @ 4a9f4b18
MODULE diaprod
! Requires key_iom_put
# if defined key_xios
!!======================================================================
!! *** MODULE diaprod ***
!! Ocean diagnostics : write ocean product diagnostics
!!=====================================================================
!! History : 3.4 ! 2012 (D. Storkey) Original code
!! 4.0 ! 2019 (D. Storkey)
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! dia_prod : calculate and write out product diagnostics
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
USE domvvl ! for thickness weighted diagnostics if key_vvl
USE eosbn2 ! equation of state (eos call)
USE phycst ! physical constants
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE in_out_manager ! I/O manager
USE iom
USE ioipsl
USE lib_mpp ! MPP library
USE timing ! preformance summary
IMPLICIT NONE
PRIVATE
PUBLIC dia_prod ! routines called by step.F90
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OPA 3.4 , NEMO Consortium (2012)
!! $Id$
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dia_prod( kt, Kmm )
!!---------------------------------------------------------------------
!! *** ROUTINE dia_prod ***
!!
!! ** Purpose : Write out product diagnostics (uT, vS etc.)
!!
!! ** Method : use iom_put
!! Product diagnostics are not thickness-weighted in this routine.
!! They should be thickness-weighted using XIOS if key_vvl is set.
!!----------------------------------------------------------------------
!!
INTEGER, INTENT( in ) :: kt ! ocean time-step index
INTEGER, INTENT( in ) :: Kmm ! ocean time level index
!!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zztmp, zztmpx, zztmpy !
!!
REAL(wp), POINTER, DIMENSION(:,:) :: z2d ! 2D workspace
REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d ! 3D workspace
REAL(wp), POINTER, DIMENSION(:,:,:) :: zrhop ! potential density
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dia_prod')
!
ALLOCATE( z2d(jpi,jpj), z3d(jpi,jpj,jpk), zrhop(jpi,jpj,jpk) )
!
IF( iom_use("urhop") .OR. iom_use("vrhop") .OR. iom_use("wrhop") &
#if ! defined key_diaar5
& .OR. iom_use("rhop") &
#endif
& ) THEN
CALL eos( ts(:,:,:,:,Kmm), z3d, zrhop ) ! now in situ and potential density
zrhop(:,:,:) = zrhop(:,:,:)-1000.e0 ! reference potential density to 1000 to avoid precision issues in rhop2 calculation
zrhop(:,:,jpk) = 0._wp
#if ! defined key_diaar5
CALL iom_put( 'rhop', zrhop )
#else
! If key_diaar5 set then there is already an iom_put call to output rhop.
! Really should be a standard diagnostics option?
#endif
ENDIF
IF( iom_use("ut") ) THEN
z3d(:,:,:) = 0.e0
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = uu(ji,jj,jk,Kmm) * 0.5 * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
END_3D
CALL iom_put( "ut", z3d ) ! product of temperature and zonal velocity at U points
ENDIF
IF( iom_use("vt") ) THEN
z3d(:,:,:) = 0.e0
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = vv(ji,jj,jk,Kmm) * 0.5 * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
END_3D
CALL iom_put( "vt", z3d ) ! product of temperature and meridional velocity at V points
ENDIF
IF( iom_use("wt") ) THEN
z3d(:,:,:) = 0.e0
DO_2D( 0, 0, 0, 0 )
z3d(ji,jj,1) = ww(ji,jj,1) * ts(ji,jj,1,jp_tem,Kmm)
END_2D
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ww(ji,jj,jk) * 0.5 * ( ts(ji,jj,jk-1,jp_tem,Kmm) + ts(ji,jj,jk,jp_tem,Kmm) )
END_3D
CALL iom_put( "wt", z3d ) ! product of temperature and vertical velocity at W points
ENDIF
IF( iom_use("us") ) THEN
z3d(:,:,:) = 0.e0
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = uu(ji,jj,jk,Kmm) * 0.5 * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
END_3D
CALL iom_put( "us", z3d ) ! product of salinity and zonal velocity at U points
ENDIF
IF( iom_use("vs") ) THEN
z3d(:,:,:) = 0.e0
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = vv(ji,jj,jk,Kmm) * 0.5 * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
END_3D
CALL iom_put( "vs", z3d ) ! product of salinity and meridional velocity at V points
ENDIF
IF( iom_use("ws") ) THEN
z3d(:,:,:) = 0.e0
DO_2D( 0, 0, 0, 0 )
z3d(ji,jj,1) = ww(ji,jj,1) * ts(ji,jj,1,jp_sal,Kmm)
END_2D
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ww(ji,jj,jk) * 0.5 * ( ts(ji,jj,jk-1,jp_sal,Kmm) + ts(ji,jj,jk,jp_sal,Kmm) )
END_3D
CALL iom_put( "ws", z3d ) ! product of salinity and vertical velocity at W points
ENDIF
IF( iom_use("uv") ) THEN
z3d(:,:,:) = 0.e0
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = 0.25 * ( uu(ji-1,jj,jk,Kmm) + uu(ji,jj,jk,Kmm) ) * ( vv(ji,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm) )
END_3D
CALL iom_put( "uv", z3d ) ! product of zonal velocity and meridional velocity at T points
ENDIF
IF( iom_use("uw") ) THEN
z3d(:,:,:) = 0.e0
DO_2D( 0, 0, 0, 0 )
z3d(ji,jj,1) = 0.5 * ( ww(ji,jj,1) + ww(ji+1,jj,1) ) * uu(ji,jj,1,Kmm)
END_2D
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = 0.25 * ( ww(ji,jj,jk) + ww(ji+1,jj,jk) ) * ( uu(ji,jj,jk-1,Kmm) + uu(ji,jj,jk,Kmm) )
END_3D
CALL iom_put( "uw", z3d ) ! product of zonal velocity and vertical velocity at UW points
ENDIF
IF( iom_use("vw") ) THEN
z3d(:,:,:) = 0.e0
DO_2D( 0, 0, 0, 0 )
z3d(ji,jj,1) = 0.5 * ( ww(ji,jj,1) + ww(ji,jj+1,1) ) * vv(ji,jj,1,Kmm)
END_2D
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = 0.25 * ( ww(ji,jj,jk) + ww(ji,jj+1,jk) ) * ( vv(ji,jj,jk-1,Kmm) + vv(ji,jj,jk,Kmm) )
END_3D
CALL iom_put( "vw", z3d ) ! product of meriodional velocity and vertical velocity at VW points
ENDIF
IF( iom_use("urhop") ) THEN
z3d(:,:,:) = 0.e0
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = uu(ji,jj,jk,Kmm) * 0.5 * ( zrhop(ji,jj,jk) + zrhop(ji+1,jj,jk) )
END_3D
CALL iom_put( "urhop", z3d ) ! product of density and zonal velocity at U points
ENDIF
IF( iom_use("vrhop") ) THEN
z3d(:,:,:) = 0.e0
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = vv(ji,jj,jk,Kmm) * 0.5 * ( zrhop(ji,jj,jk) + zrhop(ji,jj+1,jk) )
END_3D
CALL iom_put( "vrhop", z3d ) ! product of density and meridional velocity at V points
ENDIF
IF( iom_use("wrhop") ) THEN
z3d(:,:,:) = 0.e0
DO_2D( 0, 0, 0, 0 )
z3d(ji,jj,1) = ww(ji,jj,1) * zrhop(ji,jj,1)
END_2D
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ww(ji,jj,jk) * 0.5 * ( zrhop(ji,jj,jk-1) + zrhop(ji,jj,jk) )
END_3D
CALL iom_put( "wrhop", z3d ) ! product of density and vertical velocity at W points
ENDIF
!
DEALLOCATE( z2d, z3d, zrhop )
!
IF( ln_timing ) CALL timing_stop('dia_prod')
!
END SUBROUTINE dia_prod
#else
!!----------------------------------------------------------------------
!! Default option : NO diaprod
!!----------------------------------------------------------------------
LOGICAL, PUBLIC, PARAMETER :: lk_diaprod = .FALSE. ! coupled flag
CONTAINS
SUBROUTINE dia_prod( kt , Kmm) ! Empty routine
INTEGER :: kt, Kmm
!WRITE(*,*) 'dia_prod: You should not have seen this print! error?', kt
END SUBROUTINE dia_prod
#endif
!!======================================================================
END MODULE diaprod
src/OCE/DIA/diawri.F90
View file @ 4a9f4b18
......@@ -46,6 +46,7 @@ MODULE diawri
USE zdf_oce ! ocean vertical physics
USE zdfdrg ! ocean vertical physics: top/bottom friction
USE zdfmxl ! mixed layer
USE zdftke , ONLY: htau
USE zdfosm ! mixed layer
!
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
......@@ -53,6 +54,7 @@ MODULE diawri
USE dia25h ! 25h Mean output
USE iom !
USE ioipsl !
USE eosbn2
#if defined key_si3
USE ice
......@@ -124,9 +126,33 @@ CONTAINS
REAL(wp):: ze3
REAL(wp), DIMENSION(A2D( 0)) :: z2d ! 2D workspace
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z3d ! 3D workspace
CHARACTER(len=4),SAVE :: ttype , stype ! temperature and salinity type
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dia_wri')
!
IF( kt == nit000 ) THEN
IF( ln_TEOS10 ) THEN
IF ( iom_use("toce_pot") .OR. iom_use("soce_pra") .OR. iom_use("sst_pot") .OR. iom_use("sss_pra") &
& .OR. iom_use("sbt_pot") .OR. iom_use("sbs_pra") .OR. iom_use("sstgrad_pot") .OR. iom_use("sstgrad2_pot") &
& .OR. iom_use("tosmint_pot") .OR. iom_use("somint_pra")) THEN
CALL ctl_stop( 'diawri: potential temperature and practical salinity not available with ln_TEOS10' )
ELSE
ttype='con' ; stype='abs' ! teos-10 using conservative temperature and absolute salinity
ENDIF
ELSE IF ( ln_SEOS) THEN
ttype='seos' ; stype='seos' ! seos using Simplified Equation of state
ELSE
IF ( iom_use("toce_con") .OR. iom_use("soce_abs") .OR. iom_use("sst_con") .OR. iom_use("sss_abs") &
& .OR. iom_use("sbt_con") .OR. iom_use("sbs_abs") .OR. iom_use("sstgrad_con") .OR. iom_use("sstgrad2_con") &
& .OR. iom_use("tosmint_con") .OR. iom_use("somint_abs")) THEN
CALL ctl_stop( 'diawri: conservative temperature and absolute salinity not available with ln_EOS80' )
ELSE
ttype='pot' ; stype='pra' ! eos-80 using potential temperature and practical salinity
ENDIF
ENDIF
ENDIF
!
! Output the initial state and forcings
IF( ninist == 1 ) THEN
......@@ -207,25 +233,25 @@ CONTAINS
#endif
! --- tracers T&S --- !
CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) ) ! 3D temperature
CALL iom_put( "sst", ts(:,:,1,jp_tem,Kmm) ) ! surface temperature
CALL iom_put( "toce_"//ttype, ts(:,:,:,jp_tem,Kmm) ) ! 3D temperature
CALL iom_put( "sst_"//ttype, ts(:,:,1,jp_tem,Kmm) ) ! surface temperature
IF ( iom_use("sbt") ) THEN
IF ( iom_use("sbt_"//ttype) ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbkt(ji,jj)
z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
END_2D
CALL iom_put( "sbt", z2d ) ! bottom temperature
CALL iom_put( "sbt_"//ttype, z2d ) ! bottom temperature
ENDIF
CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) ) ! 3D salinity
CALL iom_put( "sss", ts(:,:,1,jp_sal,Kmm) ) ! surface salinity
IF ( iom_use("sbs") ) THEN
CALL iom_put( "soce_"//stype, ts(:,:,:,jp_sal,Kmm) ) ! 3D salinity
CALL iom_put( "sss_"//stype, ts(:,:,1,jp_sal,Kmm) ) ! surface salinity
IF ( iom_use("sbs_"//stype) ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbkt(ji,jj)
z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
END_2D
CALL iom_put( "sbs", z2d ) ! bottom salinity
CALL iom_put( "sbs_"//stype, z2d ) ! bottom salinity
ENDIF
IF( .NOT.lk_SWE ) CALL iom_put( "rhop", rhop(:,:,:) ) ! 3D potential density (sigma0)
......@@ -295,6 +321,7 @@ CONTAINS
CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef.
CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef.
CALL iom_put( "htau" , htau ) ! htau scaling
IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
......@@ -316,7 +343,7 @@ CONTAINS
ENDIF
ENDIF
IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
IF ( iom_use("sstgrad_"//ttype) .OR. iom_use("sstgrad2_"//ttype) ) THEN
DO_2D( 0, 0, 0, 0 ) ! sst gradient
zztmp = ts(ji,jj,1,jp_tem,Kmm)
zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
......@@ -324,12 +351,12 @@ CONTAINS
z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
END_2D
CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
IF ( iom_use("sstgrad") ) THEN
CALL iom_put( "sstgrad2_"//ttype, z2d ) ! square of module of sst gradient
IF ( iom_use("sstgrad_"//ttype) ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = SQRT( z2d(ji,jj) )
END_2D
CALL iom_put( "sstgrad", z2d ) ! module of sst gradient
CALL iom_put( "sstgrad_"//ttype, z2d ) ! module of sst gradient
ENDIF
ENDIF
......@@ -456,19 +483,28 @@ CONTAINS
ENDIF
IF( iom_use("tosmint") ) THEN
IF( (.NOT.l_ldfeiv_time) .AND. ( iom_use('RossRad') .OR. iom_use('RossRadlim') &
& .OR. iom_use('Tclinic_recip') .OR. iom_use('RR_GS') &
& .OR. iom_use('aeiu_2d') .OR. iom_use('aeiv_2d') ) ) THEN
CALL ldf_eiv(kt, 75.0, z2d, z3d(:,:,1), Kmm)
CALL iom_put('aeiu_2d', z2d)
CALL iom_put('aeiv_2d', z3d(:,:,1))
ENDIF
IF( iom_use("tosmint_"//ttype) ) THEN
z2d(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
END_3D
CALL iom_put( "tosmint", z2d ) ! Vertical integral of temperature
CALL iom_put( "tosmint_"//ttype, z2d ) ! Vertical integral of temperature
ENDIF
IF( iom_use("somint") ) THEN
IF( iom_use("somint_"//stype) ) THEN
z2d(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
END_3D
CALL iom_put( "somint", z2d ) ! Vertical integral of salinity
CALL iom_put( "somint_"//stype, z2d ) ! Vertical integral of salinity
ENDIF
CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
......
src/OCE/DOM/domain.F90
View file @ 4a9f4b18
......@@ -270,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
......@@ -377,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
......
src/OCE/DOM/dommsk.F90
View file @ 4a9f4b18
......@@ -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)
......@@ -207,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
......@@ -224,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
......
src/OCE/DOM/domvvl.F90
View file @ 4a9f4b18
......@@ -692,7 +692,7 @@ 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'
!
......
src/OCE/DYN/dynadv.F90
View file @ 4a9f4b18
......@@ -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)
......@@ -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_up3.F90 0 → 100644
View file @ 4a9f4b18
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/dynspg_ts.F90
View file @ 4a9f4b18
......@@ -167,7 +167,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e
REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points
REAL(wp), DIMENSION(jpi,jpj) :: zhV! fluxes
REAL(dp), DIMENSION(jpi,jpj) :: zhU! fluxes
REAL(wp), DIMENSION(jpi,jpj) :: zhU! fluxes
!!st#if defined key_qco
!!st REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v
!!st#endif
......@@ -270,6 +270,11 @@ CONTAINS
zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
END_2D
ELSE
! Ensure zhU and zhV are initialised to SOMETHING at all points to avoid referencing
! uninitialsed values in halos later on!
zhU(:,:) = 0._wp
zhV(:,:) = 0._wp
ENDIF
!
! != Add bottom stress contribution from baroclinic velocities =!
......@@ -502,13 +507,13 @@ CONTAINS
IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
!
! ! resulting flux at mid-step (not over the full domain)
DO_2D( 1, 0, 1, 1 ) ! not jpi-column
zhU(ji,jj) = e2u(ji,jj) * ua_e(ji,jj) * zhup2_e(ji,jj)
END_2D
DO_2D( 1, 1, 1, 0 ) ! not jpj-row
zhV(ji,jj) = e1v(ji,jj) * va_e(ji,jj) * zhvp2_e(ji,jj)
END_2D
!
#if defined key_agrif
! Set fluxes during predictor step to ensure volume conservation
IF( ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV )
......@@ -519,6 +524,12 @@ CONTAINS
CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V
!
ENDIF
! It seems safest to do this here since zhU and zhV are not initially calculated in halos
! by this code or by wad_Umsk, but halo values (ji-1 and jj-1) ARE required in the zhdiv
! sea level calculation. The current trunk (Feb 2024) has resolved all these issues by rewriting.
CALL lbc_lnk( 'dynspg_ts', zhU, 'U', -1._wp)
CALL lbc_lnk( 'dynspg_ts', zhV, 'V', -1._wp)
!
!
! Compute Sea Level at step jit+1
......@@ -529,15 +540,16 @@ CONTAINS
zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)
ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( ssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)
END_2D
!
CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._dp, zhU, 'U', -1._dp)
CALL lbc_lnk( 'dynspg_ts', zhV, 'V', -1._wp )
CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._dp)
!
! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
IF( ln_bdy ) CALL bdy_ssh( ssha_e )
#if defined key_agrif
CALL agrif_ssh_ts( jn )
#endif
!
! ! Sum over sub-time-steps to compute advective velocities
za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5
......@@ -556,6 +568,7 @@ CONTAINS
!
! Sea Surface Height at u-,v-points (vvl case only)
IF( .NOT.ln_linssh ) THEN
#if defined key_qcoTest_FluxForm
! ! 'key_qcoTest_FluxForm' : simple ssh average
DO_2D( 1, 0, 1, 1 )
......@@ -623,6 +636,7 @@ CONTAINS
!-- u = m+1 | h * u + delta_t' * \ h * (1-r)*g * grad_x( ssh') - h * f * k vect u + h * frc / | --!
!-- h \ / --!
!------------------------------------------------------------------------------------------------------------------------!
IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
DO_2D( 0, 0, 0, 0 )
ua_e(ji,jj) = ( un_e(ji,jj) &
......@@ -676,7 +690,7 @@ CONTAINS
ENDIF
IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
DO_2D( 0, 0, 0, 0 )
DO_2D( 0, 0, 0, 0 )
hu_e (ji,jj) = hu_0(ji,jj) + zsshu_a(ji,jj)
hur_e(ji,jj) = ssumask(ji,jj) / ( hu_e(ji,jj) + 1._wp - ssumask(ji,jj) )
hv_e (ji,jj) = hv_0(ji,jj) + zsshv_a(ji,jj)
......@@ -686,10 +700,10 @@ CONTAINS
!
IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp &
& , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp &
& , hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp )
, hu_e , 'U', 1._wp, hv_e , 'V', 1._wp &
, hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp )
ELSE
CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
ENDIF
! ! open boundaries
IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
......@@ -728,7 +742,6 @@ CONTAINS
! ! Sum sea level
pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
! ! ==================== !
END DO ! end loop !
! ! ==================== !
! -----------------------------------------------------------------------------
......@@ -1085,7 +1098,13 @@ CONTAINS
! ! Allocate time-splitting arrays
IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
!
! ! read restart when needed
! RSRH. I've just copied the following from the current NEMO main. It shouldnt affect evolution (!) but
! does help with debugging and testing!
! init some arrays for debug sette
ssha_e(:,:) = 0._wp
!
! !: restart/initialise
CALL ts_rst( nit000, 'READ' )
!
END SUBROUTINE dyn_spg_ts_init
......@@ -1165,8 +1184,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER :: ji ,jj ! dummy loop indices
REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhV
REAL(dp), DIMENSION(jpi,jpj), INTENT(in ) :: zhU
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhU, zhV
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
!!----------------------------------------------------------------------
SELECT CASE( nvor_scheme )
......@@ -1272,8 +1290,7 @@ CONTAINS
!! ** Action : ptmsk : wetting & drying t-mask
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phV, pu, pv! ocean velocities and transports
REAL(dp), DIMENSION(jpi,jpj), INTENT(inout) :: phU! ocean velocities and transports
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
!
INTEGER :: ji, jj ! dummy loop indices
......
src/OCE/DYN/dynvor.F90
View file @ 4a9f4b18
......@@ -1001,7 +1001,7 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity'
nrvm = np_RVO ! relative vorticity
ntot = np_CRV ! relative + planetary vorticity
CASE( np_FLX_c2 , np_FLX_ubs )
CASE( np_FLX_c2 , np_FLX_ubs , np_FLX_up3 )
IF(lwp) WRITE(numout,*) ' ==>>> flux form dynamics : total vorticity = Coriolis + metric term'
nrvm = np_MET ! metric term
ntot = np_CME ! Coriolis + metric term
......
src/OCE/ICB/icbrst.F90
View file @ 4a9f4b18
......@@ -24,6 +24,7 @@ MODULE icbrst
USE lib_mpp ! NEMO MPI library, lk_mpp in particular
USE netcdf ! netcdf routines for IO
USE iom
USE ioipsl, ONLY : ju2ymds ! for calendar
USE icb_oce ! define iceberg arrays
USE icbutl ! iceberg utility routines
......@@ -190,9 +191,12 @@ CONTAINS
INTEGER :: jn ! dummy loop index
INTEGER :: idg ! number of digits
INTEGER :: ix_dim, iy_dim, ik_dim, in_dim
INTEGER :: iyear, imonth, iday
REAL (wp) :: zsec
REAL (wp) :: zfjulday
CHARACTER(len=256) :: cl_path
CHARACTER(len=256) :: cl_filename
CHARACTER(len=8 ) :: cl_kt
CHARACTER(LEN=20) :: cl_kt ! ocean time-step deine as a character
CHARACTER(LEN=12 ) :: clfmt ! writing format
TYPE(iceberg), POINTER :: this
TYPE(point) , POINTER :: pt
......@@ -212,8 +216,17 @@ CONTAINS
IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/'
!
! file name
WRITE(cl_kt, '(i8.8)') kt
cl_filename = TRIM(cexper)//"_"//cl_kt//"_"//TRIM(cn_icbrst_out)
IF ( ln_rstdate ) THEN
zfjulday = fjulday + rdt / rday
IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday ) zfjulday = REAL(NINT(zfjulday),wp) ! avoid truncation error
CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )
WRITE(cl_kt, '(i4.4,2i2.2)') iyear, imonth, iday
ELSE
IF( kt > 999999999 ) THEN ; WRITE(cl_kt, * ) kt
ELSE ; WRITE(cl_kt, '(i8.8)') kt
ENDIF
ENDIF
cl_filename = TRIM(cexper)//"_"//TRIM(cl_kt)//"_"//TRIM(cn_icbrst_out)
IF( lk_mpp ) THEN
idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 ) ! how many digits to we need to write? min=4, max=9
WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg ! '(a,a,ix.x,a)'
......
src/OCE/IOM/in_out_manager.F90
View file @ 4a9f4b18
......@@ -27,6 +27,7 @@ MODULE in_out_manager
CHARACTER(lc) :: cn_ocerst_outdir !: restart output directory
LOGICAL :: ln_rstart !: start from (F) rest or (T) a restart file
LOGICAL :: ln_rst_list !: output restarts at list of times (T) or by frequency (F)
LOGICAL :: ln_rst_eos !: check equation of state used for the restart is consistent with model
INTEGER :: nn_rstctl !: control of the time step (0, 1 or 2)
INTEGER :: nn_rstssh = 0 !: hand made initilization of ssh or not (1/0)
INTEGER :: nn_it000 !: index of the first time step
......@@ -39,6 +40,7 @@ MODULE in_out_manager
INTEGER :: nn_stock !: restart file frequency
INTEGER, DIMENSION(10) :: nn_stocklist !: restart dump times
LOGICAL :: ln_mskland !: mask land points in NetCDF outputs (costly: + ~15%)
LOGICAL :: ln_rstdate !: T=> stamp output restart files with date instead of timestep
LOGICAL :: ln_cfmeta !: output additional data to netCDF files required for compliance with the CF metadata standard
LOGICAL :: ln_clobber !: clobber (overwrite) an existing file
INTEGER :: nn_chunksz !: chunksize (bytes) for NetCDF file (works only with iom_nf90 routines)
......
src/OCE/IOM/restart.F90
View file @ 4a9f4b18
......@@ -39,6 +39,7 @@ MODULE restart
!
USE in_out_manager ! I/O manager
USE iom ! I/O module
USE ioipsl, ONLY : ju2ymds ! for calendar
USE lib_mpp ! distribued memory computing library
IMPLICIT NONE
......@@ -72,6 +73,9 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time-step
!!
INTEGER :: iyear, imonth, iday
REAL (wp) :: zsec
REAL (wp) :: zfjulday
CHARACTER(LEN=20) :: clkt ! ocean time-step deine as a character
CHARACTER(LEN=50) :: clname ! ocean output restart file name
CHARACTER(lc) :: clpath ! full path to ocean output restart file
......@@ -103,8 +107,16 @@ CONTAINS
IF( kt == nitrst - 1 .OR. nn_stock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN
IF( nitrst <= nitend .AND. nitrst > 0 ) THEN
! beware of the format used to write kt (default is i8.8, that should be large enough...)
IF( nitrst > 999999999 ) THEN ; WRITE(clkt, * ) nitrst
ELSE ; WRITE(clkt, '(i8.8)') nitrst
IF ( ln_rstdate ) THEN
zfjulday = fjulday + rdt / rday
IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday ) zfjulday = REAL(NINT(zfjulday),wp) ! avoid truncation error
CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )
WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
ELSE
IF( nitrst > 999999999 ) THEN ; WRITE(clkt, * ) nitrst
ELSE ; WRITE(clkt, '(i8.8)') nitrst
ENDIF
ENDIF
! create the file
clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_ocerst_out)
......@@ -183,6 +195,9 @@ CONTAINS
#endif
ENDIF
CALL iom_rstput( kt, nitrst, numrow, 'neos' , REAL(neos)) ! equation of state
!CALL iom_rstput( kt, nitrst, numrow, 'neos' , neos , ktype = jp_i1, ldxios = lwxios) ! equation of state
IF( ln_diurnal ) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst )
IF( kt == nitrst ) THEN
IF( .NOT.lwxios ) THEN
......@@ -267,6 +282,8 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices
INTEGER :: jk
INTEGER :: id1
REAL(wp) :: zeos
REAL(wp), DIMENSION(jpi, jpj, jpk) :: w3d
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zgdept ! 3D workspace for QCO
!!----------------------------------------------------------------------
......@@ -299,6 +316,18 @@ CONTAINS
CALL iom_get( numror, jpdom_auto, 'tn', ts(:,:,:,jp_tem,Kmm) )
CALL iom_get( numror, jpdom_auto, 'sn', ts(:,:,:,jp_sal,Kmm) )
!
IF ( ln_rst_eos ) THEN
! Check equation of state used is consistent with the restart
IF( iom_varid( numror, 'neos') == -1 ) THEN
CALL ctl_stop( 'restart, rst_read: variable neos not found. STOP check that the equations of state in the restart file and in the namelist nameos are consistent and use ln_rst_eos=F')
ELSE
CALL iom_get( numror, 'neos' , zeos )
IF ( INT(zeos) /= neos ) CALL ctl_stop( 'restart, rst_read: equation of state used in restart file differs from namelist nameos')
ENDIF
ENDIF
IF( l_1st_euler ) THEN !* Euler restart (MLF only)
IF(lwp) WRITE(numout,*) ' Kbb u, v and T-S fields set to Kmm values'
uu(:,:,: ,Kbb) = uu(:,:,: ,Kmm) ! all before fields set to now values
......
src/OCE/LBC/mppini.F90
View file @ 4a9f4b18
......@@ -333,8 +333,8 @@ CONTAINS
nimpp = iimppt(ii,ij)
njmpp = ijmppt(ii,ij)
!
CALL init_doloop ! set start/end indices of do-loop, depending on the halo width value (nn_hls)
CALL init_locglo ! define now functions needed to convert indices from/to global to/from local domains
CALL init_doloop ! set start/end indices of do-loop, depending on the halo width value (nn_hls)
!
IF(lwp) THEN
WRITE(numout,*)
......@@ -1409,9 +1409,46 @@ ENDIF
Njs0 = 1+nn_hls
Nie0 = jpi-nn_hls
Nje0 = jpj-nn_hls
!
Ni_0 = Nie0 - Nis0 + 1
Nj_0 = Nje0 - Njs0 + 1
! Start and end indexes for actual coupling fields on extended grid
Nis0_ext = Nis0
Nie0_ext = Nie0
Njs0_ext = Njs0
Nje0_ext = Nje0
IF (mig(1) == 1 ) THEN
! Drag start column 1 place to the left/west
Nis0_ext=Nis0_ext-1
ENDIF
IF (mig(jpi) == jpiglo ) THEN
! Drag end column 1 place to the right/east
Nie0_ext=Nie0_ext+1
ENDIF
! RSRH we don't adjust anything in the N-S (j) direction
! since the content of the N-fold is catared for by populating values
! using lbc_lnk rather than relying on getting them from the coupler.
! This is rather confused by the fact that jpjglo has a value BIGGER
! than it did at pre 4.2... e.g. for ORCA1 it's set to 333 instead of 332
! which is rather baffling and confuses some dimensioning calculations
! which previously worked fine pre 4.2.
! Set up dimensions for old style coupling exchanges on extended grid
Ni_0_ext = Ni_0
Nj_0_ext = Nj_0
IF (mig(1) == 1 .OR. mig(jpi) == jpiglo ) THEN
! We're at the extreme left or right edge of the grid so need to cater
! for an extra column
Ni_0_ext = Ni_0_ext + 1
ENDIF
!
jpkm1 = jpk-1 ! " "
!
......
src/OCE/LDF/ldftra.F90
View file @ 4a9f4b18
......@@ -71,6 +71,8 @@ MODULE ldftra
INTEGER , PUBLIC :: nn_aei_ijk_t !: choice of time/space variation of the eiv coeff.
REAL(wp), PUBLIC :: rn_Ue !: lateral diffusive velocity [m/s]
REAL(wp), PUBLIC :: rn_Le !: lateral diffusive length [m]
INTEGER, PUBLIC :: nn_ldfeiv_shape !: shape of bounding coefficient (Treguier et al formulation only)
! ! Flag to control the type of lateral diffusive operator
INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in specification of lateral diffusion
......@@ -495,7 +497,8 @@ CONTAINS
REAL(wp) :: zah_max, zUfac ! - scalar
!!
NAMELIST/namtra_eiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv)
& nn_aei_ijk_t, rn_Ue, rn_Le ! eiv coefficient
& nn_aei_ijk_t, rn_Ue, rn_Le, & ! eiv coefficient
& nn_ldfeiv_shape
!!----------------------------------------------------------------------
!
IF(lwp) THEN ! control print
......@@ -588,6 +591,7 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( latitude, longitude, time )'
IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )'
IF(lwp) WRITE(numout,*) ' maximum allowed value: aei0 = ', aei0, ' m2/s'
IF(lwp) WRITE(numout,*) ' shape of bounding coefficient : ',nn_ldfeiv_shape
!
l_ldfeiv_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90
!
......@@ -636,14 +640,21 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: paeiu, paeiv ! eiv coefficient [m2/s]
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zfw, ze3w, zn2, z1_f20, zzaei ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zRo, zaeiw ! 2D workspace
REAL(wp) :: zfw, ze3w, zn2, z1_f20, zzaei, z2_3 ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zRo, zRo_lim, zTclinic_recip, zaeiw, zratio ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmodslp ! 3D workspace
!!----------------------------------------------------------------------
!
zn (:,:) = 0._wp ! Local initialization
zmodslp(:,:,:) = 0._wp
zhw(:,:) = 5._wp
zah(:,:) = 0._wp
zRo(:,:) = 0._wp
zRo_lim(:,:) = 0._wp
zTclinic_recip(:,:) = 0._wp
zratio(:,:) = 0._wp
zaeiw(:,:) = 0._wp
! ! Compute lateral diffusive coefficient at T-point
IF( ln_traldf_triad ) THEN
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
......@@ -670,8 +681,9 @@ CONTAINS
! eddies using the isopycnal slopes calculated in ldfslp.F :
! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w
zmodslp(ji,jj,jk) = wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
zah(ji,jj) = zah(ji,jj) + zn2 * zmodslp(ji,jj,jk) * ze3w
zhw(ji,jj) = zhw(ji,jj) + ze3w
END_3D
ENDIF
......@@ -679,17 +691,69 @@ CONTAINS
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 )
! Rossby radius at w-point taken betwenn 2 km and 40km
zRo(ji,jj) = MAX( 2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) )
zRo(ji,jj) = .4 * zn(ji,jj) / zfw
zRo_lim(ji,jj) = MAX( 2.e3 , MIN( zRo(ji,jj), 40.e3 ) )
! Compute aeiw by multiplying Ro^2 and T^-1
zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
zTclinic_recip(ji,jj) = SQRT( MAX(zah(ji,jj),0._wp) / zhw(ji,jj) ) * tmask(ji,jj,1)
zaeiw(ji,jj) = zRo_lim(ji,jj) * zRo_lim(ji,jj) * zTclinic_recip(ji,jj)
END_2D
IF( iom_use('N_2d') ) CALL iom_put('N_2d',zn(:,:)/zhw(:,:))
IF( iom_use('modslp') ) CALL iom_put('modslp',SQRT(zmodslp(:,:,:)) )
CALL iom_put('RossRad',zRo)
CALL iom_put('RossRadlim',zRo_lim)
CALL iom_put('Tclinic_recip',zTclinic_recip)
! !== Bound on eiv coeff. ==!
z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) )
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease
zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0
END_2D
z2_3 = 2._wp/3._wp
SELECT CASE(nn_ldfeiv_shape)
CASE(0) !! Standard shape applied - decrease in tropics and cap.
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease
zaeiw(ji,jj) = MIN( zzaei, paei0 )
END_2D
CASE(1) !! Abrupt cut-off on Rossby radius:
! JD : modifications here to introduce scaling by local rossby radius of deformation vs local grid scale
! arbitrary decision that GM is de-activated if local rossy radius larger than 2 times local grid scale
! based on Hallberg (2013)
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
IF ( zRo(ji,jj) >= ( 2._wp * MIN( e1t(ji,jj), e2t(ji,jj) ) ) ) THEN
! TODO : use a version of zRo that integrates over a few time steps ?
zaeiw(ji,jj) = 0._wp
ELSE
zaeiw(ji,jj) = MIN( zaeiw(ji,jj), paei0 )
ENDIF
END_2D
CASE(2) !! Rossby radius ramp type 1:
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zratio(ji,jj) = zRo(ji,jj)/MIN(e1t(ji,jj),e2t(ji,jj))
zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, z2_3*(2._wp - zratio(ji,jj)) ) ) * paei0 )
END_2D
CALL iom_put('RR_GS',zratio)
CASE(3) !! Rossby radius ramp type 2:
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zratio(ji,jj) = MIN(e1t(ji,jj),e2t(ji,jj))/zRo(ji,jj)
zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, z2_3*( zratio(ji,jj) - 0.5_wp ) ) ) * paei0 )
END_2D
CASE(4) !! bathymetry ramp:
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, 0.001*(ht_0(ji,jj) - 2000._wp) ) ) * paei0 )
END_2D
CASE(5) !! Rossby radius ramp type 1 applied to Treguier et al coefficient rather than cap:
!! Note the ramp is RR/GS=[2.0,1.0] (not [2.0,0.5] as for cases 2,3) and we ramp up
!! to 5% of the Treguier et al coefficient, aiming for peak values of around 100m2/s
!! at high latitudes rather than 2000m2/s which is what you get in eORCA025 with an
!! uncapped coefficient.
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zratio(ji,jj) = zRo(ji,jj)/MIN(e1t(ji,jj),e2t(ji,jj))
zaeiw(ji,jj) = MAX( 0._wp, MIN( 1._wp, 2._wp - zratio(ji,jj) ) ) * 0.05 * zaeiw(ji,jj)
zaeiw(ji,jj) = MIN( zaeiw(ji,jj), paei0 )
END_2D
CALL iom_put('RR_GS',zratio)
CASE DEFAULT
CALL ctl_stop('ldf_eiv: Unrecognised option for nn_ldfeiv_shape.')
END SELECT
IF( nn_hls == 1 ) CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1.0_wp ) ! lateral boundary condition
!
DO_2D( 0, 0, 0, 0 )
......
src/OCE/SBC/cpl_oasis3.F90
View file @ 4a9f4b18
......@@ -30,20 +30,35 @@ MODULE cpl_oasis3
USE xios ! I/O server
#endif
USE par_oce ! ocean parameters
USE cpl_rnf_1d, ONLY: nn_cpl_river ! Variables used in 1D river outflow
USE dom_oce ! ocean space and time domain
USE in_out_manager ! I/O manager
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE lib_mpp
#if defined key_agrif || ! defined key_mpi_off
USE MPI
#endif
IMPLICIT NONE
PRIVATE
#if ! defined key_oasis3
! Dummy interface to oasis_get if not using oasis
INTERFACE oasis_get
MODULE PROCEDURE oasis_get_1d, oasis_get_2d
END INTERFACE
#endif
PUBLIC cpl_init
PUBLIC cpl_define
PUBLIC cpl_snd
PUBLIC cpl_rcv
PUBLIC cpl_rcv_1d
PUBLIC cpl_freq
PUBLIC cpl_finalize
INTEGER, PARAMETER :: localRoot = 0
INTEGER, PUBLIC :: OASIS_Rcv = 1 !: return code if received field
INTEGER, PUBLIC :: OASIS_idle = 0 !: return code if nothing done by oasis
INTEGER :: ncomp_id ! id returned by oasis_init_comp
......@@ -67,9 +82,13 @@ MODULE cpl_oasis3
INTEGER :: nrcv ! total number of fields received
INTEGER :: nsnd ! total number of fields sent
INTEGER :: ncplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
INTEGER, PUBLIC, PARAMETER :: nmaxfld=62 ! Maximum number of coupling fields
INTEGER, PUBLIC, PARAMETER :: nmaxfld=100 ! Maximum number of coupling fields
INTEGER, PUBLIC, PARAMETER :: nmaxcat=5 ! Maximum number of coupling fields
INTEGER, PUBLIC, PARAMETER :: nmaxcpl=5 ! Maximum number of coupling fields
LOGICAL, PARAMETER :: ltmp_wapatch = .FALSE. ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define
LOGICAL, PARAMETER :: ltmp_landproc = .TRUE. ! patch to restrict coupled area to non halo cells
TYPE, PUBLIC :: FLD_CPL !: Type for coupling field information
LOGICAL :: laction ! To be coupled or not
......@@ -79,11 +98,16 @@ MODULE cpl_oasis3
INTEGER, DIMENSION(nmaxcat,nmaxcpl) :: nid ! Id of the field (no more than 9 categories and 9 extrena models)
INTEGER :: nct ! Number of categories in field
INTEGER :: ncplmodel ! Maximum number of models to/from which this variable may be sent/received
INTEGER :: dimensions ! Number of dimensions of coupling field
END TYPE FLD_CPL
TYPE(FLD_CPL), DIMENSION(nmaxfld), PUBLIC :: srcv, ssnd !: Coupling fields
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: exfld ! Temporary buffer for receiving
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: exfld_ext ! Temporary buffer for receiving with wrap points
INTEGER :: ishape_ext(4) ! shape of 2D arrays passed to PSMILe extended for wrap points in weights data
INTEGER :: ishape(4) ! shape of 2D arrays passed to PSMILe
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -103,6 +127,7 @@ CONTAINS
!!--------------------------------------------------------------------
CHARACTER(len = *), INTENT(in ) :: cd_modname ! model name as set in namcouple file
INTEGER , INTENT( out) :: kl_comm ! local communicator of the model
INTEGER :: error
!!--------------------------------------------------------------------
! WARNING: No write in numout in this routine
......@@ -123,6 +148,7 @@ CONTAINS
IF( nerror /= OASIS_Ok ) &
CALL oasis_abort (ncomp_id, 'cpl_init','Failure in oasis_get_localcomm' )
!
END SUBROUTINE cpl_init
......@@ -138,13 +164,26 @@ CONTAINS
INTEGER, INTENT(in) :: krcv, ksnd ! Number of received and sent coupling fields
INTEGER, INTENT(in) :: kcplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
!
INTEGER :: id_part
INTEGER :: id_part_0d ! Partition for 0d fields
INTEGER :: id_part_rnf_1d ! Partition for 1d river outflow fields
INTEGER :: id_part_2d ! Partition for 2d fields
INTEGER :: id_part_2d_ext ! Partition for 2d fields extended for old (pre vn4.2) style remapping weights!
INTEGER :: id_part_temp ! Temperary partition used to choose either 0d or 1d partitions
INTEGER :: paral(5) ! OASIS3 box partition
INTEGER :: ishape(4) ! shape of arrays passed to PSMILe
INTEGER :: paral_ext(5) ! OASIS3 box partition extended
INTEGER :: ishape0d1d(2) ! Shape of 0D or 1D arrays passed to PSMILe.
INTEGER :: var_nodims(2) ! Number of coupling field dimensions.
! var_nodims(1) is redundant from OASIS3-MCT vn4.0 onwards
! but retained for backward compatibility.
! var_nodims(2) is the number of fields in a bundle
! or 1 for unbundled fields (bundles are not yet catered for
! in NEMO hence we default to 1).
INTEGER :: ji,jc,jm ! local loop indicees
LOGICAL :: llenddef ! should we call xios_oasis_enddef and oasis_enddef?
CHARACTER(LEN=64) :: zclname
CHARACTER(LEN=2) :: cli2
INTEGER :: i_offset ! Used in calculating offset for extended partition.
!!--------------------------------------------------------------------
IF(lwp) WRITE(numout,*)
......@@ -173,10 +212,13 @@ CONTAINS
ishape(2) = Ni_0
ishape(3) = 1
ishape(4) = Nj_0
!
ishape0d1d(1) = 0
ishape0d1d(2) = 0 !
! ... Allocate memory for data exchange
!
ALLOCATE(exfld(Ni_0, Nj_0), stat = nerror) ! allocate only inner domain (without halos)
ALLOCATE(exfld(Ni_0, Nj_0), stat = nerror) ! allocate full domain (without halos)
IF( nerror > 0 ) THEN
CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld') ; RETURN
ENDIF
......@@ -198,7 +240,96 @@ CONTAINS
WRITE(numout,*) ' multiexchg: Njs0, Nje0, njmpp =', Njs0, Nje0, njmpp
ENDIF
CALL oasis_def_partition ( id_part, paral, nerror, Ni0glo*Nj0glo ) ! global number of points, excluding halos
! We still set up the new vn4.2 style box partition for reference, though it doesn't actually get used,
! we can easily swap back to it if we ever manage to successfully generate vn4.2 compatible weights, or introduce
! RTL controls to distinguish between onl and new style weights.
CALL oasis_def_partition ( id_part_2d, paral, nerror, Ni0glo*Nj0glo ) ! global number of points, excluding halos
! RSRH Set up 2D box partition for compatibility with existing weights files
! so we don't have to generate and manage multiple sets of weights purely because of
! the changes to nemo 4.2+ code!
! This is just a hack for global cyclic models for the time being
Ni0glo_ext = jpiglo
Nj0glo_ext = Nj0glo +1 ! We can't use jpjglo here because for some reason at 4.2 this is bigger
! than at 4.0.... e.g. for ORCA1 it is 333 when it should only be 332!
! RSRH extended shapes for old style dimensioning. Allows backwards compatibility with existing weights files,
! which the new code DOES NOT, causing headaches not only for users but also for management of weights files.
ishape_ext(:) = ishape(:)
IF (mig(1) == 1 .OR. mig(jpi)==jpiglo) THEN
! Extra columns in PEs dealing with wrap points
ishape_ext(2) = ishape_ext(2) + 1
ENDIF
! Workout any extra offset in the i dimension
IF (mig(1) == 1 ) THEN
i_offset = mig0(nn_hls)
ELSE
i_offset = mig(nn_hls)
ENDIF
ALLOCATE(exfld_ext(ishape_ext(2), ishape_ext(4)), stat = nerror) ! allocate full domain (with wrap pts)
IF( nerror > 0 ) THEN
CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld_ext') ; RETURN
ENDIF
! Now we have the appropriate dimensions, we can set up the partition array for the old-style extended grid
paral_ext(1) = 2 ! box partitioning
paral_ext(2) = (Ni0glo_ext * mjg0(nn_hls)) + i_offset ! NEMO lower left corner global offset, with wrap pts
paral_ext(3) = Ni_0_ext ! local extent in i, including halos
paral_ext(4) = Nj_0_ext ! local extent in j, including halos
paral_ext(5) = Ni0glo_ext ! global extent in x, including halos
IF( sn_cfctl%l_oasout ) THEN
WRITE(numout,*) ' multiexchg: paral_ext (1:5)', paral_ext, jpiglo, jpjglo, Ni0glo_ext, Nj0glo_ext
WRITE(numout,*) ' multiexchg: Ni_0_ext, Nj_0_ext i_offset =', Ni_0_ext, Nj_0_ext, i_offset
WRITE(numout,*) ' multiexchg: Nis0_ext, Nie0_ext =', Nis0_ext, Nie0_ext
WRITE(numout,*) ' multiexchg: Njs0_ext, Nje0_ext =', Njs0_ext, Nje0_ext
ENDIF
! Define our extended grid
CALL oasis_def_partition ( id_part_2d_ext, paral_ext, nerror, Ni0glo_ext*Nj0glo_ext )
! OK so now we should have a usable 2d partition for fields defined WITH redundant points.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! A special partition is needed for 0D fields
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
paral(1) = 0 ! serial partitioning
paral(2) = 0
IF ( mpprank == 0) THEN
paral(3) = 1 ! Size of array to couple (scalar)
ELSE
paral(3) = 0 ! Dummy size for PE's not involved
END IF
paral(4) = 0
paral(5) = 0
CALL oasis_def_partition ( id_part_0d, paral, nerror, 1 )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Another special partition is needed for 1D river routing fields
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
paral(1) = 0 ! serial partitioning
paral(2) = 0
IF ( mpprank == 0) THEN
paral(3) = nn_cpl_river ! Size of array to couple (vector)
ELSE
paral(3) = 0 ! Dummy size for PE's not involved
END IF
paral(4) = 0
paral(5) = 0
CALL oasis_def_partition ( id_part_rnf_1d, paral, nerror, nn_cpl_river )
!
! ... Announce send variables.
!
......@@ -232,14 +363,24 @@ CONTAINS
ENDIF
#endif
IF( sn_cfctl%l_oasout ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_Out
CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part , (/ 2, 1 /), &
& OASIS_Out , ishape , OASIS_REAL, nerror )
flush(numout)
IF( sn_cfctl%l_oasout ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_Out
CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part_2d_ext , (/ 2, 1 /), &
& OASIS_Out , ishape_ext , OASIS_REAL, nerror )
IF( nerror /= OASIS_Ok ) THEN
WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname)
CALL oasis_abort ( ssnd(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' )
ENDIF
IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "variable defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "variable NOT defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "put variable defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "put variable NOT defined in the namcouple"
END DO
END DO
ENDIF
......@@ -276,15 +417,57 @@ CONTAINS
zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname)
ENDIF
#endif
IF( sn_cfctl%l_oasout ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_In
CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part , (/ 2, 1 /), &
& OASIS_In , ishape , OASIS_REAL, nerror )
flush(numout)
! Define 0D (Greenland or Antarctic ice mass) or 1D (river outflow) coupling fields
IF (srcv(ji)%dimensions <= 1) THEN
var_nodims(1) = 1
var_nodims(2) = 1 ! Modify this value to cater for bundled fields.
IF (mpprank == 0) THEN
IF (srcv(ji)%dimensions == 0) THEN
! If 0D then set temporary variables to 0D components
id_part_temp = id_part_0d
ishape0d1d(2) = 1
ELSE
! If 1D then set temporary variables to river outflow components
id_part_temp = id_part_rnf_1d
ishape0d1d(2)= nn_cpl_river
END IF
CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part_temp , var_nodims, &
OASIS_In , ishape0d1d(1:2) , OASIS_REAL, nerror )
ELSE
! Dummy call to keep OASIS3-MCT happy.
CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part_0d , var_nodims, &
OASIS_In , ishape0d1d(1:2) , OASIS_REAL, nerror )
END IF
ELSE
! It's a "normal" 2D (or pseudo 3D) coupling field.
! ... Set the field dimension and bundle count
var_nodims(1) = 2
var_nodims(2) = 1 ! Modify this value to cater for bundled fields.
CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part_2d_ext , var_nodims, &
OASIS_In , ishape_ext , OASIS_REAL, nerror )
ENDIF
IF( nerror /= OASIS_Ok ) THEN
WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname)
CALL oasis_abort ( srcv(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' )
ENDIF
IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "variable defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "variable NOT defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "get variable defined in the namcouple"
IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "get variable NOT defined in the namcouple"
END DO
END DO
......@@ -330,13 +513,17 @@ CONTAINS
INTEGER :: jc,jm ! local loop index
!!--------------------------------------------------------------------
!
! snd data to OASIS3
!
DO jc = 1, ssnd(kid)%nct
DO jm = 1, ssnd(kid)%ncplmodel
IF( ssnd(kid)%nid(jc,jm) /= -1 ) THEN ! exclude halos from data sent to oasis
CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(Nis0:Nie0, Njs0:Nje0,jc), kinfo )
! The field is "put" directly, using appropriate start/end indexing - i.e. we don't
! copy it to an intermediate buffer.
CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc), kinfo )
IF ( sn_cfctl%l_oasout ) THEN
IF ( kinfo == OASIS_Sent .OR. kinfo == OASIS_ToRest .OR. &
......@@ -346,10 +533,11 @@ CONTAINS
WRITE(numout,*) 'oasis_put: ivarid ', ssnd(kid)%nid(jc,jm)
WRITE(numout,*) 'oasis_put: kstep ', kstep
WRITE(numout,*) 'oasis_put: info ', kinfo
WRITE(numout,*) ' - Minimum value is ', MINVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Maximum value is ', MAXVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Sum value is ', SUM(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Minimum value is ', MINVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) ' - Maximum value is ', MAXVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) ' - Sum value is ', SUM(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) '****************'
CALL FLUSH(numout)
ENDIF
ENDIF
......@@ -357,6 +545,7 @@ CONTAINS
ENDDO
ENDDO
!
END SUBROUTINE cpl_snd
......@@ -375,13 +564,16 @@ CONTAINS
INTEGER , INTENT( out) :: kinfo ! OASIS3 info argument
!!
INTEGER :: jc,jm ! local loop index
LOGICAL :: llaction, ll_1st
LOGICAL :: llaction, ll_1st, lrcv
!!--------------------------------------------------------------------
!
! receive local data from OASIS3 on every process
!
kinfo = OASIS_idle
lrcv=.FALSE.
!
DO jc = 1, srcv(kid)%nct
ll_1st = .TRUE.
......@@ -389,7 +581,7 @@ CONTAINS
IF( srcv(kid)%nid(jc,jm) /= -1 ) THEN
CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld, kinfo )
CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld_ext, kinfo )
llaction = kinfo == OASIS_Recvd .OR. kinfo == OASIS_FromRest .OR. &
& kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut
......@@ -398,14 +590,18 @@ CONTAINS
& WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
IF( llaction ) THEN ! data received from oasis do not include halos
! but DO still cater for wrap columns when using pre vn4.2 compatible remapping weights.
kinfo = OASIS_Rcv
lrcv=.TRUE.
IF( ll_1st ) THEN
pdata(Nis0:Nie0,Njs0:Nje0,jc) = exfld(:,:) * pmask(Nis0:Nie0,Njs0:Nje0,jm)
pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc) = exfld_ext(:,:) * pmask(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jm)
ll_1st = .FALSE.
ELSE
pdata(Nis0:Nie0,Njs0:Nje0,jc) = pdata(Nis0:Nie0,Njs0:Nje0,jc) &
& + exfld(:,:) * pmask(Nis0:Nie0,Njs0:Nje0,jm)
pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc) = pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc) &
& + exfld_ext(:,:) * pmask(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jm)
ENDIF
IF ( sn_cfctl%l_oasout ) THEN
......@@ -414,10 +610,11 @@ CONTAINS
WRITE(numout,*) 'oasis_get: ivarid ' , srcv(kid)%nid(jc,jm)
WRITE(numout,*) 'oasis_get: kstep', kstep
WRITE(numout,*) 'oasis_get: info ', kinfo
WRITE(numout,*) ' - Minimum value is ', MINVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Maximum value is ', MAXVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Sum value is ', SUM(pdata(Nis0:Nie0,Njs0:Nje0,jc))
WRITE(numout,*) ' - Minimum value is ', MINVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) ' - Maximum value is ', MAXVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) ' - Sum value is ', SUM(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
WRITE(numout,*) '****************'
CALL FLUSH(numout)
ENDIF
ENDIF
......@@ -426,15 +623,130 @@ CONTAINS
ENDDO
!--- we must call lbc_lnk to fill the halos that where not received.
IF( .NOT. ll_1st ) THEN
CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )
ENDIF
ENDDO
! RSRH I've changed this since:
! 1) it seems multi cat fields may not be properly updated in halos when called on a per
! category basis(?)
! 2) it's more efficient to have a single call (and simpler to understand) to update ALL
! categories at the same time!
!--- Call lbc_lnk to populate halos of received fields.
IF (lrcv) then
CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,:), srcv(kid)%clgrid, srcv(kid)%nsgn )
endif
!
END SUBROUTINE cpl_rcv
SUBROUTINE cpl_rcv_1d( kid, kstep, pdata, nitems, kinfo )
!!---------------------------------------------------------------------
!! *** ROUTINE cpl_rcv_1d ***
!!
!! ** Purpose : - A special version of cpl_rcv to deal exclusively with
!! receipt of 0D or 1D fields.
!! The fields are recieved into a 1D array buffer which is simply a
!! dynamically sized sized array (which may be of size 1)
!! of 0 dimensional fields. This allows us to pass miltiple 0D
!! fields via a single put/get operation.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: nitems ! Number of 0D items to recieve
! during this get operation. i.e.
! The size of the 1D array in which
! 0D items are passed.
INTEGER , INTENT(in ) :: kid ! ID index of the incoming
! data.
INTEGER , INTENT(in ) :: kstep ! ocean time-step in seconds
REAL(wp), INTENT(inout) :: pdata(1:nitems) ! The original value(s),
! unchanged if nothing is
! received
INTEGER , INTENT( out) :: kinfo ! OASIS3 info argument
!!
REAL(wp) :: recvfld(1:nitems) ! Local receive field buffer
INTEGER :: jc,jm ! local loop index
INTEGER :: ierr
LOGICAL :: llaction
INTEGER :: MPI_WORKING_PRECISION
INTEGER :: number_to_print
!!--------------------------------------------------------------------
!
! receive local data from OASIS3 on every process
!
kinfo = OASIS_idle
!
! 0D and 1D fields won't have categories or any other form of "pseudo level"
! so we only cater for a single set of values and thus don't bother
! with a loop over the jc index
jc = 1
DO jm = 1, srcv(kid)%ncplmodel
IF( srcv(kid)%nid(jc,jm) /= -1 ) THEN
IF ( ( srcv(kid)%dimensions <= 1) .AND. (mpprank == 0) ) THEN
! Since there is no concept of data decomposition for zero
! dimension fields, they must only be exchanged through the master PE,
! unlike "normal" 2D field cases where every PE is involved.
CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, recvfld, kinfo )
llaction = kinfo == OASIS_Recvd .OR. kinfo == OASIS_FromRest .OR. &
kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut
IF ( sn_cfctl%l_oasout ) WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , &
llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
IF ( llaction ) THEN
kinfo = OASIS_Rcv
pdata(1:nitems) = recvfld(1:nitems)
IF( sn_cfctl%l_oasout ) THEN
number_to_print = 10
IF ( nitems < number_to_print ) number_to_print = nitems
WRITE(numout,*) '****************'
WRITE(numout,*) 'oasis_get: Incoming ', srcv(kid)%clname
WRITE(numout,*) 'oasis_get: ivarid ' , srcv(kid)%nid(jc,jm)
WRITE(numout,*) 'oasis_get: kstep', kstep
WRITE(numout,*) 'oasis_get: info ', kinfo
WRITE(numout,*) ' - Minimum Value is ', MINVAL(pdata(:))
WRITE(numout,*) ' - Maximum value is ', MAXVAL(pdata(:))
WRITE(numout,*) ' - Start of data is ', pdata(1:number_to_print)
WRITE(numout,*) '****************'
CALL FLUSH(numout)
ENDIF
ENDIF
ENDIF
ENDIF
ENDDO
#if defined key_mpi_off
CALL oasis_abort( ncomp_id, "cpl_rcv_1d", "Unable to use mpi_bcast with key_mpi_off in your list of NEMO keys." )
#else
! Set the precision that we want to broadcast using MPI_BCAST
SELECT CASE( wp )
CASE( sp )
MPI_WORKING_PRECISION = MPI_REAL ! Single precision
CASE( dp )
MPI_WORKING_PRECISION = MPI_DOUBLE_PRECISION ! Double precision
CASE default
CALL oasis_abort( ncomp_id, "cpl_rcv_1d", "Could not find precision for coupling 0d or 1d field" )
END SELECT
! We have to broadcast (potentially) received values from PE 0 to all
! the others. If no new data has been received we're just
! broadcasting the existing values but there's no more efficient way
! to deal with that w/o NEMO adopting a UM-style test mechanism
! to determine active put/get timesteps.
CALL mpi_bcast( pdata, nitems, MPI_WORKING_PRECISION, localRoot, mpi_comm_oce, ierr )
#endif
!
END SUBROUTINE cpl_rcv_1d
INTEGER FUNCTION cpl_freq( cdfieldname )
!!---------------------------------------------------------------------
......@@ -500,6 +812,7 @@ CONTAINS
!!----------------------------------------------------------------------
!
DEALLOCATE( exfld )
DEALLOCATE( exfld_ext )
IF(nstop == 0) THEN
CALL oasis_terminate( nerror )
ELSE
......@@ -543,7 +856,7 @@ CONTAINS
SUBROUTINE oasis_def_var(k1,cd1,k2,k3,k4,k5,k6,k7)
CHARACTER(*), INTENT(in ) :: cd1
INTEGER , INTENT(in ) :: k2,k3(2),k4,k5(2,2),k6
INTEGER , INTENT(in ) :: k2,k3(2),k4,k5(*),k6
INTEGER , INTENT( out) :: k1,k7
k1 = -1 ; k7 = -1
WRITE(numout,*) 'oasis_def_var: Error you sould not be there...', cd1
......@@ -563,13 +876,21 @@ CONTAINS
WRITE(numout,*) 'oasis_put: Error you sould not be there...'
END SUBROUTINE oasis_put
SUBROUTINE oasis_get(k1,k2,p1,k3)
SUBROUTINE oasis_get_2d(k1,k2,p1,k3)
REAL(wp), DIMENSION(:,:), INTENT( out) :: p1
INTEGER , INTENT(in ) :: k1,k2
INTEGER , INTENT( out) :: k3
p1(1,1) = -1. ; k3 = -1
WRITE(numout,*) 'oasis_get: Error you sould not be there...'
END SUBROUTINE oasis_get
WRITE(numout,*) 'oasis_get_2d: Error you sould not be there...'
END SUBROUTINE oasis_get_2d
SUBROUTINE oasis_get_1d(k1,k2,p1,k3)
REAL(wp), DIMENSION(:) , INTENT( out) :: p1
INTEGER , INTENT(in ) :: k1,k2
INTEGER , INTENT( out) :: k3
p1(1) = -1. ; k3 = -1
WRITE(numout,*) 'oasis_get_1d: Error you sould not be there...'
END SUBROUTINE oasis_get_1d
SUBROUTINE oasis_get_freqs(k1,k5,k2,k3,k4)
INTEGER , INTENT(in ) :: k1,k2
......
src/OCE/SBC/cpl_rnf_1d.F90 0 → 100644
View file @ 4a9f4b18
MODULE cpl_rnf_1d
!!======================================================================
!! *** MODULE cpl_rnf_1d ***
!! Ocean forcing: River runoff passed from the atmosphere using
!! 1D array. One value per river.
!!=====================================================================
!! History : ?.? ! 2018-01 (D. Copsey) Initial setup
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! cpl_rnf_1d_init : runoffs initialisation
!!----------------------------------------------------------------------
#if defined key_oasis3
USE mod_oasis ! OASIS3-MCT module
#endif
USE timing ! Timing
USE in_out_manager ! I/O units
USE lib_mpp ! MPP library
USE iom
USE dom_oce ! Domain sizes (for grid box area e1e2t)
USE sbc_oce ! Surface boundary condition: ocean fields
USE lib_fortran, ONLY: DDPDD
IMPLICIT NONE
PRIVATE
PUBLIC cpl_rnf_1d_init ! routine called in nemo_init
PUBLIC cpl_rnf_1d_to_2d ! routine called in sbccpl.F90
TYPE, PUBLIC :: RIVERS_DATA !: Storage for river outflow data
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: river_number !: River outflow number
REAL(wp), ALLOCATABLE, DIMENSION(:) :: river_area ! 1D array listing areas of each river outflow (m2)
COMPLEX(wp), ALLOCATABLE, DIMENSION(:) :: river_area_c ! Comlex version of river_area for use in bit reproducible sums (m2)
END TYPE RIVERS_DATA
! DB: see https://forge.ipsl.jussieu.fr/nemo/ticket/1865
#if defined key_agrif
TYPE(RIVERS_DATA), PUBLIC :: rivers !: River data
#else
TYPE(RIVERS_DATA), PUBLIC, TARGET :: rivers !: River data
#endif
INTEGER, PUBLIC :: nn_cpl_river ! Maximum number of rivers being passed through the coupler
INTEGER, PUBLIC :: runoff_id ! OASIS coupling id used in oasis_get command
LOGICAL :: ln_print_river_info ! Diagnostic prints of river coupling information
CONTAINS
SUBROUTINE cpl_rnf_1d_init
!!----------------------------------------------------------------------
!! *** SUBROUTINE cpl_rnf_1d_init ***
!!
!! ** Purpose : - Read in file for river outflow numbers.
!! Calculate 2D area of river outflow points.
!! Called from nemo_init (nemogcm.F90).
!!
!!----------------------------------------------------------------------
!! namelist variables
!!-------------------
CHARACTER(len=200) :: file_riv_number !: Filename for river numbers
INTEGER :: ios ! Local integer output status for namelist read
INTEGER :: inum
INTEGER :: ii, jj, rr !: Loop indices
INTEGER :: max_river
REAL(wp), DIMENSION(jpi,jpj) :: river_number ! 2D array containing the river outflow numbers
NAMELIST/nam_cpl_rnf_1d/file_riv_number, nn_cpl_river, ln_print_river_info
!!----------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('cpl_rnf_1d_init')
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'cpl_rnf_1d_init : initialization of river runoff coupling'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
IF(lwp) CALL flush(numout)
!! REWIND(numnam_cfg)
! Read the namelist
READ ( numnam_ref, nam_cpl_rnf_1d, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_cpl_rnf_1d in reference namelist' )
READ ( numnam_cfg, nam_cpl_rnf_1d, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_cpl_rnf_1d in configuration namelist' )
IF(lwm) WRITE ( numond, nam_cpl_rnf_1d )
! ! Parameter control and print
IF(lwp) WRITE(numout,*) ' '
IF(lwp) WRITE(numout,*) ' Namelist nam_cpl_rnf_1d : Coupled runoff using 1D array'
IF(lwp) WRITE(numout,*) ' Input file that contains river numbers = ',file_riv_number
IF(lwp) WRITE(numout,*) ' Maximum number of rivers to couple = ',nn_cpl_river
IF(lwp) WRITE(numout,*) ' Print river information = ',ln_print_river_info
IF(lwp) WRITE(numout,*) ' '
IF(lwp) CALL flush(numout)
! Assign space for river numbers
ALLOCATE( rivers%river_number( jpi, jpj ) )
! Read the river numbers from netcdf file
CALL iom_open (file_riv_number , inum )
CALL iom_get ( inum, jpdom_global, 'river_number', river_number, cd_type = 'T', psgn = 1._wp)
CALL iom_close( inum )
! Convert from a real array to an integer array
max_river=0
DO ii = 1, jpi
DO jj = 1, jpj
rivers%river_number(ii,jj) = INT(river_number(ii,jj))
IF ( rivers%river_number(ii,jj) > max_river ) THEN
max_river = rivers%river_number(ii,jj)
END IF
END DO
END DO
! Print out the largest river number
IF ( ln_print_river_info .AND. lwp) THEN
WRITE(numout,*) 'Maximum river number in input file = ',max_river
CALL flush(numout)
END IF
! Get the area of each river outflow
ALLOCATE( rivers%river_area( nn_cpl_river ) )
ALLOCATE( rivers%river_area_c( nn_cpl_river ) )
rivers%river_area_c(:) = CMPLX( 0.e0, 0.e0, wp )
DO ii = Nis0, Nie0
DO jj = Njs0, Nje0
IF ( tmask_i(ii,jj) > 0.5 ) THEN ! This makes sure we are not at a duplicated point (at north fold or east-west cyclic point)
IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river ) THEN
! Add the area of each grid box (e1e2t) into river_area_c using DDPDD which should maintain bit reproducibility (needs to be checked)
CALL DDPDD( CMPLX( e1e2t(ii,jj), 0.e0, wp ), rivers%river_area_c(rivers%river_number(ii,jj) ) )
END IF
END IF
END DO
END DO
! Use mpp_sum to add together river areas on other processors
CALL mpp_sum( 'cpl_rnf_1d', rivers%river_area_c )
! Convert from complex number to real
rivers%river_area(:) = REAL(rivers%river_area_c(:),wp)
IF ( ln_print_river_info .AND. lwp) THEN
WRITE(numout,*) 'Area of rivers 1 to 10 are ',rivers%river_area(1:10)
WRITE(numout,*) 'Maximum river area = ',MAXVAL(rivers%river_area)
WRITE(numout,*) 'Minimum river area = ',MINVAL(rivers%river_area)
WRITE(numout,*) 'Sum of river areas = ',SUM(rivers%river_area)
CALL flush(numout)
END IF
IF ( MINVAL(rivers%river_area) <= 0 ) THEN
WRITE(numout,*) 'ERROR: There is at least one river with a river outflow area of zero. Please check your file_riv_number file'
WRITE(numout,*) 'that all the allocated river numbers are at ocean points as defined by the bathymetry file with no river'
WRITE(numout,*) 'numbers within the north fold or wraparound points.'
DO rr = 1,nn_cpl_river
IF ( rivers%river_area(rr) <= 0 ) THEN
WRITE(numout,*) 'The first river with this problem is river number ',rr
CALL ctl_stop ( 'STOP', 'ERROR: There is at least one river with a river outflow area of zero. See stdout.')
END IF
END DO
END IF
END SUBROUTINE cpl_rnf_1d_init
SUBROUTINE cpl_rnf_1d_to_2d( runoff_1d )
!!----------------------------------------------------------------------
!! *** SUBROUTINE cpl_rnf_1d_to_2d ***
!!
!! ** Purpose : - Convert river outflow from 1D array (passed from the
!! atmosphere) to the 2D NEMO runoff field.
!! Called from sbc_cpl_ice_flx (sbccpl.F90).
!!
!!----------------------------------------------------------------------
REAL , INTENT(in ) :: runoff_1d(nn_cpl_river) ! River runoff. One value per river.
INTEGER :: ii, jj, rr ! Loop indices
LOGICAL :: found_nan
! Convert the 1D total runoff per river to 2D runoff flux by
! dividing by the area of each runoff zone.
DO ii = 1, jpi
DO jj = 1, jpj
IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river .AND. tmask_i(ii,jj) > 0.5 ) THEN
rnf(ii,jj) = runoff_1d(rivers%river_number(ii,jj)) / rivers%river_area(rivers%river_number(ii,jj))
ELSE
rnf(ii,jj) = 0.0
END IF
END DO
END DO
IF ( ln_print_river_info ) THEN
WRITE(numout,*) 'SUM of river runoff on 1D array = ',SUM(runoff_1d)
WRITE(numout,*) 'SUM of river runoff on 2D array = ',SUM(rnf)
found_nan = .false.
DO rr = 1, nn_cpl_river
IF (rr <= 10) THEN
WRITE(numout,*) 'River number ',rr,' has total runoff=',runoff_1d(rr),' and area=',rivers%river_area(rr),' making runoff flux=',runoff_1d(rr)/rivers%river_area(rr)
END IF
IF (runoff_1d(rr) /= runoff_1d(rr)) THEN
IF (.NOT. found_nan) THEN
WRITE(numout,*) 'WARNING: NAN found at river number ',rr
found_nan = .true.
END IF
END IF
END DO
END IF
END SUBROUTINE cpl_rnf_1d_to_2d
END MODULE cpl_rnf_1d