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
  • nemo/nemo
  • sparonuz/nemo
  • hatfield/nemo
  • extdevs/nemo
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1355 additions and 907 deletions
src/OCE/SBC/cpl_oasis3.F90
View file @ ee6ce8c5
......@@ -419,6 +419,7 @@ CONTAINS
IF( .NOT. ll_1st ) THEN
CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )
ENDIF
!!clem: mettre T instead of clgrid
ENDDO
!
......
src/OCE/SBC/sbc_ice.F90
View file @ ee6ce8c5
......@@ -50,8 +50,8 @@ MODULE sbc_ice
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qcn_ice !: heat conduction flux in the layer below surface [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qtr_ice_top !: solar flux transmitted below the ice surface [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_ice !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtau_ice !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_ice !: atmos-ice u-stress. T-pts [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtau_ice !: atmos-ice v-stress. T-pts [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp_ice !: sublimation - precip over sea ice [kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: topmelt !: category topmelt
......
src/OCE/SBC/sbc_oce.F90
View file @ ee6ce8c5
......@@ -103,34 +103,36 @@ MODULE sbc_oce
INTEGER , PUBLIC :: ncpl_qsr_freq = 0 !: qsr coupling frequency per days from atmosphere (used by top)
!
!! !! now ! before !!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau , utau_b !: sea surface i-stress (ocean referential) [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtau , vtau_b !: sea surface j-stress (ocean referential) [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_icb, vtau_icb !: sea surface (i,j)-stress used by icebergs [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: taum !: module of sea surface stress (at T-point) [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau !: sea surface i-stress (ocean referential) T-pt [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtau !: sea surface j-stress (ocean referential) T-pt [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utauU , utau_b !: sea surface i-stress (ocean referential) U-pt [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtauV , vtau_b !: sea surface j-stress (ocean referential) V-pt [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_icb, vtau_icb !: sea surface (i,j)-stress used by icebergs [N/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: taum !: module of sea surface stress (at T-point) [N/m2]
!! wndm is used compute surface gases exchanges in ice-free ocean or leads
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wndm !: wind speed module at T-point (=|U10m-Uoce|) [m/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhoa !: air density at "rn_zu" m above the sea [kg/m3]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qsr !: sea heat flux: solar [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qns , qns_b !: sea heat flux: non solar [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qsr_tot !: total solar heat flux (over sea and ice) [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qns_tot !: total non solar heat flux (over sea and ice) [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp , emp_b !: freshwater budget: volume flux [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx , sfx_b !: salt flux [PSS.kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp_tot !: total E-P over ocean and ice [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fmmflx !: freshwater budget: freezing/melting [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnf , rnf_b !: river runoff [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwficb , fwficb_b !: iceberg melting [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wndm !: wind speed module at T-point (=|U10m-Uoce|) [m/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhoa !: air density at "rn_zu" m above the sea [kg/m3]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qsr !: sea heat flux: solar [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qns , qns_b !: sea heat flux: non solar [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qsr_tot !: total solar heat flux (over sea and ice) [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qns_tot !: total non solar heat flux (over sea and ice) [W/m2]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp , emp_b !: freshwater budget: volume flux [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx , sfx_b !: salt flux [PSS.kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp_tot !: total E-P over ocean and ice [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fmmflx !: freshwater budget: freezing/melting [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnf , rnf_b !: river runoff [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwficb , fwficb_b !: iceberg melting [Kg/m2/s]
!!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sbc_tsc, sbc_tsc_b !: sbc content trend [K.m/s] jpi,jpj,jpts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qsr_hc , qsr_hc_b !: heat content trend due to qsr flux [K.m/s] jpi,jpj,jpk
!!
!!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tprecip !: total precipitation [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sprecip !: solid precipitation [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fr_i !: ice fraction = 1 - lead fraction (between 0 to 1)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: atm_co2 !: atmospheric pCO2 [ppm]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: cloud_fra !: cloud cover (fraction of cloud in a gridcell) [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tprecip !: total precipitation [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sprecip !: solid precipitation [Kg/m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fr_i !: ice fraction = 1 - lead fraction (between 0 to 1)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: atm_co2 !: atmospheric pCO2 [ppm]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: cloud_fra !: cloud cover (fraction of cloud in a gridcell) [-]
!!---------------------------------------------------------------------
!! ABL Vertical Domain size
......@@ -177,8 +179,8 @@ CONTAINS
!!---------------------------------------------------------------------
ierr(:) = 0
!
ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , taum(jpi,jpj) , &
& vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) )
ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , utauU(jpi,jpj) , taum(jpi,jpj) , &
& vtau(jpi,jpj) , vtau_b(jpi,jpj) , vtauV(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) )
!
ALLOCATE( qns_tot(jpi,jpj) , qns (jpi,jpj) , qns_b(jpi,jpj), &
& qsr_tot(jpi,jpj) , qsr (jpi,jpj) , &
......@@ -205,9 +207,10 @@ CONTAINS
END FUNCTION sbc_oce_alloc
!!clem => this subroutine is never used in nemo
SUBROUTINE sbc_tau2wnd
!!---------------------------------------------------------------------
!! *** ROUTINE sbc_tau2wnd ***
!! *** ROUTINE ***
!!
!! ** Purpose : Estimation of wind speed as a function of wind stress
!!
......@@ -217,17 +220,14 @@ CONTAINS
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3
REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient
REAL(wp) :: ztx, zty, ztau, zcoef ! temporary variables
REAL(wp) :: ztau, zcoef ! temporary variables
INTEGER :: ji, jj ! dummy indices
!!---------------------------------------------------------------------
zcoef = 0.5 / ( zrhoa * zcdrag )
DO_2D( 0, 0, 0, 0 )
ztx = utau(ji-1,jj ) + utau(ji,jj)
zty = vtau(ji ,jj-1) + vtau(ji,jj)
ztau = SQRT( ztx * ztx + zty * zty )
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
ztau = SQRT( utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj) )
wndm(ji,jj) = SQRT ( ztau * zcoef ) * tmask(ji,jj,1)
END_2D
CALL lbc_lnk( 'sbc_oce', wndm(:,:) , 'T', 1.0_wp )
!
END SUBROUTINE sbc_tau2wnd
......
src/OCE/SBC/sbcblk.F90
View file @ ee6ce8c5
......@@ -493,7 +493,7 @@ CONTAINS
!! the stress is assumed to be in the (i,j) mesh referential
!!
!! ** Action : defined at each time-step at the air-sea interface
!! - utau, vtau i- and j-component of the wind stress
!! - utau, vtau i- and j-component of the wind stress at T-point
!! - taum wind stress module at T-point
!! - wndm wind speed module at T-point over free ocean or leads in presence of sea-ice
!! - qns, qsr non-solar and solar heat fluxes
......@@ -611,10 +611,10 @@ CONTAINS
END SUBROUTINE sbc_blk
SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair, & ! inp
& pslp , pst , pu , pv, & ! inp
& puatm, pvatm, pdqsr , pdqlw , & ! inp
& ptsk , pssq , pcd_du, psen, plat, pevp ) ! out
SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair, & ! <<= in
& pslp , pst , pu , pv, & ! <<= in
& puatm, pvatm, pdqsr , pdqlw , & ! <<= in
& ptsk , pssq , pcd_du, psen, plat, pevp ) ! =>> out
!!---------------------------------------------------------------------
!! *** ROUTINE blk_oce_1 ***
!!
......@@ -657,7 +657,6 @@ CONTAINS
#if defined key_cyclone
REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point
#endif
REAL(wp), DIMENSION(jpi,jpj) :: ztau_i, ztau_j ! wind stress components at T-point
REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]
REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean
REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean
......@@ -695,22 +694,20 @@ CONTAINS
#else
! ... scalar wind module at T-point (not masked)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
END_2D
#endif
! ----------------------------------------------------------------------------- !
! I Solar FLUX !
! ----------------------------------------------------------------------------- !
! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave
zztmp = 1. - albo
! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle
IF( ln_dm2dc ) THEN
qsr(:,:) = zztmp * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1)
qsr(:,:) = ( 1._wp - albo ) * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1)
ELSE
qsr(:,:) = zztmp * pdqsr(:,:) * tmask(:,:,1)
qsr(:,:) = ( 1._wp - albo ) * pdqsr(:,:) * tmask(:,:,1)
ENDIF
! ----------------------------------------------------------------------------- !
! II Turbulent FLUXES !
! ----------------------------------------------------------------------------- !
......@@ -718,69 +715,62 @@ CONTAINS
! specific humidity at SST
pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) )
! Backup "bulk SST" and associated spec. hum.
IF( ln_skin_cs .OR. ln_skin_wl ) THEN
!! Backup "bulk SST" and associated spec. hum.
zztmp1(:,:) = zsspt(:,:)
zztmp2(:,:) = pssq(:,:)
zztmp2(:,:) = pssq (:,:)
ENDIF
!! Time to call the user-selected bulk parameterization for
!! == transfer coefficients ==! Cd, Ch, Ce at T-point, and more...
SELECT CASE( nblk )
! transfer coefficients (Cd, Ch, Ce at T-point, and more)
SELECT CASE( nblk ) ! user-selected bulk parameterization
!
CASE( np_NCAR )
CALL turb_ncar ( rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
& nb_iter=nn_iter_algo )
!
CASE( np_COARE_3p0 )
CALL turb_coare3p0( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_COARE_3p6 )
CALL turb_coare3p6( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_ECMWF )
CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_ANDREAS )
CALL turb_andreas ( rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo )
!
CASE DEFAULT
CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk parameterizaton selected' )
!
END SELECT
IF( iom_use('Cd_oce') ) CALL iom_put("Cd_oce", zcd_oce * tmask(:,:,1))
IF( iom_use('Ce_oce') ) CALL iom_put("Ce_oce", zce_oce * tmask(:,:,1))
IF( iom_use('Ch_oce') ) CALL iom_put("Ch_oce", zch_oce * tmask(:,:,1))
! outputs
IF( iom_use('Cd_oce') ) CALL iom_put( "Cd_oce", zcd_oce * tmask(:,:,1) )
IF( iom_use('Ce_oce') ) CALL iom_put( "Ce_oce", zce_oce * tmask(:,:,1) )
IF( iom_use('Ch_oce') ) CALL iom_put( "Ch_oce", zch_oce * tmask(:,:,1) )
!! LB: mainly here for debugging purpose:
IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ptair-rt0) * tmask(:,:,1)) ! potential temperature at z=zt
IF( iom_use('q_zt') ) CALL iom_put("q_zt", pqair * tmask(:,:,1)) ! specific humidity "
IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (theta_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu
IF( iom_use('q_zu') ) CALL iom_put("q_zu", q_zu * tmask(:,:,1)) ! specific humidity "
IF( iom_use('ssq') ) CALL iom_put("ssq", pssq * tmask(:,:,1)) ! saturation specific humidity at z=0
IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu
IF( iom_use('theta_zt') ) CALL iom_put( "theta_zt", (ptair-rt0) * tmask(:,:,1) ) ! potential temperature at z=zt
IF( iom_use('q_zt') ) CALL iom_put( "q_zt", pqair * tmask(:,:,1) ) ! specific humidity "
IF( iom_use('theta_zu') ) CALL iom_put( "theta_zu", (theta_zu -rt0) * tmask(:,:,1) ) ! potential temperature at z=zu
IF( iom_use('q_zu') ) CALL iom_put( "q_zu", q_zu * tmask(:,:,1) ) ! specific humidity "
IF( iom_use('ssq') ) CALL iom_put( "ssq", pssq * tmask(:,:,1) ) ! saturation specific humidity at z=0
IF( iom_use('wspd_blk') ) CALL iom_put( "wspd_blk", zU_zu * tmask(:,:,1) ) ! bulk wind speed at z=zu
! In the presence of sea-ice we do not use the cool-skin/warm-layer update of zsspt, pssq & ptsk from turb_*()
IF( ln_skin_cs .OR. ln_skin_wl ) THEN
!! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zsspt, pssq & ptsk:
WHERE ( fr_i(:,:) > 0.001_wp )
! sea-ice present, we forget about the update, using what we backed up before call to turb_*()
zsspt(:,:) = zztmp1(:,:)
pssq(:,:) = zztmp2(:,:)
pssq (:,:) = zztmp2(:,:)
END WHERE
! apply potential temperature increment to abolute SST
ptsk(:,:) = ptsk(:,:) + ( zsspt(:,:) - zztmp1(:,:) )
......@@ -809,10 +799,10 @@ CONTAINS
END_2D
CALL BULK_FORMULA( rn_zu, zsspt(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
& wndm(:,:), zU_zu(:,:), pslp(:,:), rhoa(:,:), &
& taum(:,:), psen(:,:), plat(:,:), &
& pEvap=pevp(:,:), pfact_evap=rn_efac )
& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
& wndm(:,:), zU_zu(:,:), pslp(:,:), rhoa(:,:), &
& taum(:,:), psen(:,:), plat(:,:), &
& pEvap=pevp(:,:), pfact_evap=rn_efac )
psen(:,:) = psen(:,:) * tmask(:,:,1)
plat(:,:) = plat(:,:) * tmask(:,:,1)
......@@ -821,57 +811,42 @@ CONTAINS
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF( wndm(ji,jj) > 0._wp ) THEN
zztmp = taum(ji,jj) / wndm(ji,jj)
zztmp = taum(ji,jj) / wndm(ji,jj)
#if defined key_cyclone
ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
utau(ji,jj) = zztmp * zwnd_i(ji,jj)
vtau(ji,jj) = zztmp * zwnd_j(ji,jj)
#else
ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
utau(ji,jj) = zztmp * pwndi(ji,jj)
vtau(ji,jj) = zztmp * pwndj(ji,jj)
#endif
ELSE
ztau_i(ji,jj) = 0._wp
ztau_j(ji,jj) = 0._wp
utau(ji,jj) = 0._wp
vtau(ji,jj) = 0._wp
ENDIF
END_2D
IF( ln_crt_fbk ) THEN ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715)
zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp ) ! set the max value of Stau corresponding to a wind of 3 m/s (<0)
DO_2D( 0, 1, 0, 1 ) ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop
zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax ) ! stau (<0) must be smaller than zstmax
ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj ) + pu(ji,jj) ) - puatm(ji,jj) )
ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) )
taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) )
DO_2D( 0, 0, 0, 0 )
zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax ) * tmask(ji,jj,1) ! stau (<0) must be smaller than zstmax
utau(ji,jj) = utau(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj ) + pu(ji,jj) ) - puatm(ji,jj) )
vtau(ji,jj) = vtau(ji,jj) + zstau * ( 0.5_wp * ( pv(ji ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) )
taum(ji,jj) = SQRT( utau(ji,jj) * utau(ji,jj) + vtau(ji,jj) * vtau(ji,jj) )
END_2D
CALL lbc_lnk( 'sbcblk', utau, 'T', -1._wp, vtau, 'T', -1._wp, taum, 'T', 1._wp )
ENDIF
! ... utau, vtau at U- and V_points, resp.
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
! Note that coastal wind stress is not used in the code... so this extra care has no effect
DO_2D( 0, 0, 0, 0 ) ! start loop at 2, in case ln_crt_fbk = T
utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj ) ) &
& * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji ,jj+1) ) &
& * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
END_2D
IF( ln_crt_fbk ) THEN
CALL lbc_lnk( 'sbcblk', utau, 'U', -1._wp, vtau, 'V', -1._wp, taum, 'T', 1._wp )
ELSE
CALL lbc_lnk( 'sbcblk', utau, 'U', -1._wp, vtau, 'V', -1._wp )
ENDIF
! Saving open-ocean wind-stress (module and components) on T-points:
CALL iom_put( "taum_oce", taum(:,:)*tmask(:,:,1) ) ! output wind stress module
!#LB: These 2 lines below mostly here for 'STATION_ASF' test-case, otherwize "utau" (U-grid) and vtau" (V-grid) does the job in: [DYN/dynatf.F90])
CALL iom_put( "utau_oce", ztau_i(:,:)*tmask(:,:,1) ) ! utau at T-points!
CALL iom_put( "vtau_oce", ztau_j(:,:)*tmask(:,:,1) ) ! vtau at T-points!
! Saving open-ocean wind-stress (module and components)
CALL iom_put( "taum_oce", taum(:,:) ) ! wind stress module
! ! LB: These 2 lines below mostly here for 'STATION_ASF' test-case
CALL iom_put( "utau_oce", utau(:,:) ) ! utau
CALL iom_put( "vtau_oce", vtau(:,:) ) ! vtau
IF(sn_cfctl%l_prtctl) THEN
CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ', mask1=tmask )
CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ', mask1=tmask )
CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, &
& tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask )
CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=tmask, &
& tab2d_2=vtau , clinfo2=' vtau : ', mask2=tmask )
CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd : ', mask1=tmask )
ENDIF
!
......@@ -896,8 +871,8 @@ CONTAINS
!! at the ocean surface at each time step knowing Cd, Ch, Ce and
!! atmospheric variables (from ABL or external data)
!!
!! ** Outputs : - utau : i-component of the stress at U-point (N/m2)
!! - vtau : j-component of the stress at V-point (N/m2)
!! ** Outputs : - utau : i-component of the stress at T-point (N/m2)
!! - vtau : j-component of the stress at T-point (N/m2)
!! - taum : Wind stress module at T-point (N/m2)
!! - wndm : Wind speed module at T-point (m/s)
!! - qsr : Solar heat flux over the ocean (W/m2)
......@@ -1015,11 +990,16 @@ CONTAINS
REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zootm_su ! sea-ice surface mean temperature
REAL(wp) :: zztmp1, zztmp2 ! temporary scalars
REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsipt ! temporary array
REAL(wp) :: zztmp ! temporary scalars
REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsipt ! temporary array
REAL(wp), DIMENSION(jpi,jpj) :: zmsk00 ! O% concentration ice mask
!!---------------------------------------------------------------------
!
! treshold for outputs
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , fr_i(ji,jj) - 1.e-6_wp ) ) ! 1 if ice, 0 if no ice
END_2D
! ------------------------------------------------------------ !
! Wind module relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
......@@ -1032,9 +1012,9 @@ CONTAINS
zsipt(:,:) = theta_exner( ptsui(:,:), pslp(:,:) )
! sea-ice <-> atmosphere bulk transfer coefficients
SELECT CASE( nblk_ice )
CASE( np_ice_cst )
SELECT CASE( nblk_ice ) ! user-selected bulk parameterization
!
CASE( np_ice_cst )
! Constant bulk transfer coefficients over sea-ice:
Cd_ice(:,:) = rn_Cd_i
Ch_ice(:,:) = rn_Ch_i
......@@ -1042,62 +1022,45 @@ CONTAINS
! no height adjustment, keeping zt values:
theta_zu_i(:,:) = ptair(:,:)
q_zu_i(:,:) = pqair(:,:)
CASE( np_ice_an05 ) ! calculate new drag from Lupkes(2015) equations
!
CASE( np_ice_an05 ) ! from Andreas(2005) equations
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
CALL turb_ice_an05( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
CASE( np_ice_lu12 )
!
CASE( np_ice_lu12 ) ! from Lupkes(2012) equations
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
CALL turb_ice_lu12( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
CASE( np_ice_lg15 ) ! calculate new drag from Lupkes(2015) equations
!
CASE( np_ice_lg15 ) ! from Lupkes and Gryanik (2015) equations
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
CALL turb_ice_lg15( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
!
END SELECT
IF( iom_use('Cd_ice').OR.iom_use('Ce_ice').OR.iom_use('Ch_ice').OR.iom_use('taum_ice').OR.iom_use('utau_ice').OR.iom_use('vtau_ice') ) &
& ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice !
IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice*ztmp)
IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice*ztmp)
IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice*ztmp)
IF( ln_blk ) THEN
! ---------------------------------------------------- !
! Wind stress relative to nonmoving ice ( U10m ) !
! ---------------------------------------------------- !
! supress moving ice in wind stress computation as we don't know how to do it properly...
DO_2D( 0, 1, 0, 1 ) ! at T point
zztmp1 = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
putaui(ji,jj) = zztmp1 * pwndi(ji,jj)
pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zztmp = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
putaui(ji,jj) = zztmp * pwndi(ji,jj)
pvtaui(ji,jj) = zztmp * pwndj(ji,jj)
END_2D
!#LB: saving the module, and x-y components, of the ai wind-stress at T-points: NOT weighted by the ice concentration !!!
IF(iom_use('taum_ice')) CALL iom_put('taum_ice', SQRT( putaui*putaui + pvtaui*pvtaui )*ztmp )
!#LB: These 2 lines below mostly here for 'STATION_ASF' test-case, otherwize "utau_oi" (U-grid) and vtau_oi" (V-grid) does the job in: [ICE/icedyn_rhg_evp.F90])
IF(iom_use('utau_ice')) CALL iom_put("utau_ice", putaui*ztmp) ! utau at T-points!
IF(iom_use('vtau_ice')) CALL iom_put("vtau_ice", pvtaui*ztmp) ! vtau at T-points!
!
DO_2D( 0, 0, 0, 0 ) ! U & V-points (same as ocean).
!#LB: QUESTION?? so SI3 expects wind stress vector to be provided at U & V points? Not at T-points ?
! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology
zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) )
zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) )
putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj ) )
pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji ,jj+1) )
END_2D
CALL lbc_lnk( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )
! outputs
! LB: not weighted by the ice concentration
IF( iom_use('taum_ice') ) CALL iom_put( 'taum_ice', SQRT( putaui*putaui + pvtaui*pvtaui ) * zmsk00 )
! LB: These 2 lines below mostly here for 'STATION_ASF' test-case
IF( iom_use('utau_ice') ) CALL iom_put( "utau_ice", putaui * zmsk00 )
IF( iom_use('vtau_ice') ) CALL iom_put( "vtau_ice", pvtaui * zmsk00 )
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ', mask1=umask &
& , tab2d_2=pvtaui , clinfo2=' pvtaui : ', mask2=vmask )
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ', mask1=tmask &
& , tab2d_2=pvtaui , clinfo2=' pvtaui : ', mask2=tmask )
ELSE ! ln_abl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
......@@ -1105,10 +1068,15 @@ CONTAINS
pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
END_2D
pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ; ! more accurate way to obtain ssq !
pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! more accurate way to obtain ssq
ENDIF ! ln_blk / ln_abl
!
! outputs
IF( iom_use('Cd_ice') ) CALL iom_put( "Cd_ice", Cd_ice * zmsk00 )
IF( iom_use('Ce_ice') ) CALL iom_put( "Ce_ice", Ce_ice * zmsk00 )
IF( iom_use('Ch_ice') ) CALL iom_put( "Ch_ice", Ch_ice * zmsk00 )
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ', mask1=tmask )
!
END SUBROUTINE blk_ice_1
......
src/OCE/SBC/sbccpl.F90
View file @ ee6ce8c5
......@@ -68,15 +68,15 @@ MODULE sbccpl
INTEGER, PARAMETER :: jpr_otx1 = 1 ! 3 atmosphere-ocean stress components on grid 1
INTEGER, PARAMETER :: jpr_oty1 = 2 !
INTEGER, PARAMETER :: jpr_otz1 = 3 !
INTEGER, PARAMETER :: jpr_otx2 = 4 ! 3 atmosphere-ocean stress components on grid 2
INTEGER, PARAMETER :: jpr_oty2 = 5 !
INTEGER, PARAMETER :: jpr_otz2 = 6 !
!!$ INTEGER, PARAMETER :: jpr_otx2 = 4 ! 3 atmosphere-ocean stress components on grid 2
!!$ INTEGER, PARAMETER :: jpr_oty2 = 5 !
!!$ INTEGER, PARAMETER :: jpr_otz2 = 6 !
INTEGER, PARAMETER :: jpr_itx1 = 7 ! 3 atmosphere-ice stress components on grid 1
INTEGER, PARAMETER :: jpr_ity1 = 8 !
INTEGER, PARAMETER :: jpr_itz1 = 9 !
INTEGER, PARAMETER :: jpr_itx2 = 10 ! 3 atmosphere-ice stress components on grid 2
INTEGER, PARAMETER :: jpr_ity2 = 11 !
INTEGER, PARAMETER :: jpr_itz2 = 12 !
!!$ INTEGER, PARAMETER :: jpr_itx2 = 10 ! 3 atmosphere-ice stress components on grid 2
!!$ INTEGER, PARAMETER :: jpr_ity2 = 11 !
!!$ INTEGER, PARAMETER :: jpr_itz2 = 12 !
INTEGER, PARAMETER :: jpr_qsroce = 13 ! Qsr above the ocean
INTEGER, PARAMETER :: jpr_qsrice = 14 ! Qsr above the ice
INTEGER, PARAMETER :: jpr_qsrmix = 15
......@@ -128,9 +128,9 @@ MODULE sbccpl
INTEGER, PARAMETER :: jpr_isf = 60
INTEGER, PARAMETER :: jpr_icb = 61
INTEGER, PARAMETER :: jpr_ts_ice = 62 ! Sea ice surface temp
!!INTEGER, PARAMETER :: jpr_qtrice = 63 ! Transmitted solar thru sea-ice
INTEGER, PARAMETER :: jpr_qtrice = 63 ! Transmitted solar thru sea-ice
INTEGER, PARAMETER :: jprcv = 62 ! total number of fields received
INTEGER, PARAMETER :: jprcv = 63 ! total number of fields received
INTEGER, PARAMETER :: jps_fice = 1 ! ice fraction sent to the atmosphere
INTEGER, PARAMETER :: jps_toce = 2 ! ocean temperature
......@@ -194,7 +194,7 @@ MODULE sbccpl
& sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr
! ! Received from the atmosphere
TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, &
& sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf, sn_rcv_ts_ice
& sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf, sn_rcv_ts_ice, sn_rcv_qtrice
TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf
! ! Send to waves
TYPE(FLD_C) :: sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev
......@@ -202,7 +202,6 @@ MODULE sbccpl
TYPE(FLD_C) :: sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, &
& sn_rcv_wstrf, sn_rcv_wdrag, sn_rcv_charn, sn_rcv_taw, sn_rcv_bhd, sn_rcv_tusd, sn_rcv_tvsd
! ! Other namelist parameters
!! TYPE(FLD_C) :: sn_rcv_qtrice
INTEGER :: nn_cplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
LOGICAL :: ln_usecplmask ! use a coupling mask file to merge data received from several models
! -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
......@@ -281,7 +280,7 @@ CONTAINS
& sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum , sn_rcv_wstrf , &
& sn_rcv_charn , sn_rcv_taw , sn_rcv_bhd , sn_rcv_tusd , sn_rcv_tvsd, &
& sn_rcv_wdrag , sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf , sn_rcv_cal , &
& sn_rcv_iceflx, sn_rcv_co2 , sn_rcv_icb , sn_rcv_isf , sn_rcv_ts_ice, & !!, sn_rcv_qtrice
& sn_rcv_iceflx, sn_rcv_co2 , sn_rcv_icb , sn_rcv_isf , sn_rcv_ts_ice, sn_rcv_qtrice, &
& sn_rcv_mslp
!!---------------------------------------------------------------------
......@@ -313,7 +312,6 @@ CONTAINS
WRITE(numout,*)' surface stress = ', TRIM(sn_rcv_tau%cldes ), ' (', TRIM(sn_rcv_tau%clcat ), ')'
WRITE(numout,*)' - referential = ', sn_rcv_tau%clvref
WRITE(numout,*)' - orientation = ', sn_rcv_tau%clvor
WRITE(numout,*)' - mesh = ', sn_rcv_tau%clvgrd
WRITE(numout,*)' non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
WRITE(numout,*)' solar heat flux = ', TRIM(sn_rcv_qsr%cldes ), ' (', TRIM(sn_rcv_qsr%clcat ), ')'
WRITE(numout,*)' non-solar heat flux = ', TRIM(sn_rcv_qns%cldes ), ' (', TRIM(sn_rcv_qns%clcat ), ')'
......@@ -323,7 +321,7 @@ CONTAINS
WRITE(numout,*)' iceberg = ', TRIM(sn_rcv_icb%cldes ), ' (', TRIM(sn_rcv_icb%clcat ), ')'
WRITE(numout,*)' ice shelf = ', TRIM(sn_rcv_isf%cldes ), ' (', TRIM(sn_rcv_isf%clcat ), ')'
WRITE(numout,*)' sea ice heat fluxes = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
!! WRITE(numout,*)' transmitted solar thru sea-ice = ', TRIM(sn_rcv_qtrice%cldes), ' (', TRIM(sn_rcv_qtrice%clcat), ')'
WRITE(numout,*)' transmitted solar thru sea-ice = ', TRIM(sn_rcv_qtrice%cldes), ' (', TRIM(sn_rcv_qtrice%clcat), ')'
WRITE(numout,*)' atm co2 = ', TRIM(sn_rcv_co2%cldes ), ' (', TRIM(sn_rcv_co2%clcat ), ')'
WRITE(numout,*)' Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')'
WRITE(numout,*)' surface waves:'
......@@ -358,7 +356,7 @@ CONTAINS
WRITE(numout,*)' - mesh = ', sn_snd_crtw%clvgrd
ENDIF
IF( lwp .AND. ln_wave) THEN ! control print
WRITE(numout,*)' surface waves:'
WRITE(numout,*)' surface waves:'
WRITE(numout,*)' Significant wave heigth = ', TRIM(sn_rcv_hsig%cldes ), ' (', TRIM(sn_rcv_hsig%clcat ), ')'
WRITE(numout,*)' Wave to oce energy flux = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'
WRITE(numout,*)' Surface Stokes drift grid u = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')'
......@@ -368,8 +366,8 @@ CONTAINS
WRITE(numout,*)' Stress frac adsorbed by waves = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')'
WRITE(numout,*)' Neutral surf drag coefficient = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')'
WRITE(numout,*)' Charnock coefficient = ', TRIM(sn_rcv_charn%cldes ), ' (', TRIM(sn_rcv_charn%clcat ), ')'
WRITE(numout,*)' Transport associated to Stokes drift grid u = ', TRIM(sn_rcv_tusd%cldes ), ' (', TRIM(sn_rcv_tusd%clcat ), ')'
WRITE(numout,*)' Transport associated to Stokes drift grid v = ', TRIM(sn_rcv_tvsd%cldes ), ' (', TRIM(sn_rcv_tvsd%clcat ), ')'
WRITE(numout,*)' Transport associated to Stokes drift u = ', TRIM(sn_rcv_tusd%cldes ), ' (', TRIM(sn_rcv_tusd%clcat ), ')'
WRITE(numout,*)' Transport associated to Stokes drift v = ', TRIM(sn_rcv_tvsd%cldes ), ' (', TRIM(sn_rcv_tvsd%clcat ), ')'
WRITE(numout,*)' Bernouilli pressure head = ', TRIM(sn_rcv_bhd%cldes ), ' (', TRIM(sn_rcv_bhd%clcat ), ')'
WRITE(numout,*)'Wave to ocean momentum flux and Net wave-supported stress = ', TRIM(sn_rcv_taw%cldes ), ' (', TRIM(sn_rcv_taw%clcat ), ')'
WRITE(numout,*)' Surface current to waves = ', TRIM(sn_snd_crtw%cldes ), ' (', TRIM(sn_snd_crtw%clcat ), ')'
......@@ -399,87 +397,26 @@ CONTAINS
srcv(jpr_otx1)%clname = 'O_OTaux1' ! 1st ocean component on grid ONE (T or U)
srcv(jpr_oty1)%clname = 'O_OTauy1' ! 2nd - - - -
srcv(jpr_otz1)%clname = 'O_OTauz1' ! 3rd - - - -
srcv(jpr_otx2)%clname = 'O_OTaux2' ! 1st ocean component on grid TWO (V)
srcv(jpr_oty2)%clname = 'O_OTauy2' ! 2nd - - - -
srcv(jpr_otz2)%clname = 'O_OTauz2' ! 3rd - - - -
!
srcv(jpr_itx1)%clname = 'O_ITaux1' ! 1st ice component on grid ONE (T, F, I or U)
srcv(jpr_ity1)%clname = 'O_ITauy1' ! 2nd - - - -
srcv(jpr_itz1)%clname = 'O_ITauz1' ! 3rd - - - -
srcv(jpr_itx2)%clname = 'O_ITaux2' ! 1st ice component on grid TWO (V)
srcv(jpr_ity2)%clname = 'O_ITauy2' ! 2nd - - - -
srcv(jpr_itz2)%clname = 'O_ITauz2' ! 3rd - - - -
!
! Vectors: change of sign at north fold ONLY if on the local grid
IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM( sn_rcv_tau%cldes ) == 'oce and ice' &
.OR. TRIM( sn_rcv_tau%cldes ) == 'mixed oce-ice' ) THEN ! avoid working with the atmospheric fields if they are not coupled
!
IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' ) srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
! ! Set grid and action
SELECT CASE( TRIM( sn_rcv_tau%clvgrd ) ) ! 'T', 'U,V', 'U,V,I', 'U,V,F', 'T,I', 'T,F', or 'T,U,V'
CASE( 'T' )
srcv(jpr_otx1:jpr_itz2)%clgrid = 'T' ! oce and ice components given at T-point
srcv(jpr_otx1:jpr_otz1)%laction = .TRUE. ! receive oce components on grid 1
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1
CASE( 'U,V' )
srcv(jpr_otx1:jpr_otz1)%clgrid = 'U' ! oce components given at U-point
srcv(jpr_otx2:jpr_otz2)%clgrid = 'V' ! and V-point
srcv(jpr_itx1:jpr_itz1)%clgrid = 'U' ! ice components given at U-point
srcv(jpr_itx2:jpr_itz2)%clgrid = 'V' ! and V-point
srcv(jpr_otx1:jpr_itz2)%laction = .TRUE. ! receive oce and ice components on both grid 1 & 2
CASE( 'U,V,T' )
srcv(jpr_otx1:jpr_otz1)%clgrid = 'U' ! oce components given at U-point
srcv(jpr_otx2:jpr_otz2)%clgrid = 'V' ! and V-point
srcv(jpr_itx1:jpr_itz1)%clgrid = 'T' ! ice components given at T-point
srcv(jpr_otx1:jpr_otz2)%laction = .TRUE. ! receive oce components on grid 1 & 2
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1 only
CASE( 'U,V,I' )
srcv(jpr_otx1:jpr_otz1)%clgrid = 'U' ! oce components given at U-point
srcv(jpr_otx2:jpr_otz2)%clgrid = 'V' ! and V-point
srcv(jpr_itx1:jpr_itz1)%clgrid = 'I' ! ice components given at I-point
srcv(jpr_otx1:jpr_otz2)%laction = .TRUE. ! receive oce components on grid 1 & 2
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1 only
CASE( 'U,V,F' )
srcv(jpr_otx1:jpr_otz1)%clgrid = 'U' ! oce components given at U-point
srcv(jpr_otx2:jpr_otz2)%clgrid = 'V' ! and V-point
srcv(jpr_itx1:jpr_itz1)%clgrid = 'F' ! ice components given at F-point
srcv(jpr_otx1:jpr_otz2)%laction = .TRUE. ! receive oce components on grid 1 & 2
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1 only
CASE( 'T,I' )
srcv(jpr_otx1:jpr_itz2)%clgrid = 'T' ! oce and ice components given at T-point
srcv(jpr_itx1:jpr_itz1)%clgrid = 'I' ! ice components given at I-point
srcv(jpr_otx1:jpr_otz1)%laction = .TRUE. ! receive oce components on grid 1
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1
CASE( 'T,F' )
srcv(jpr_otx1:jpr_itz2)%clgrid = 'T' ! oce and ice components given at T-point
srcv(jpr_itx1:jpr_itz1)%clgrid = 'F' ! ice components given at F-point
srcv(jpr_otx1:jpr_otz1)%laction = .TRUE. ! receive oce components on grid 1
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1
CASE( 'T,U,V' )
srcv(jpr_otx1:jpr_otz1)%clgrid = 'T' ! oce components given at T-point
srcv(jpr_itx1:jpr_itz1)%clgrid = 'U' ! ice components given at U-point
srcv(jpr_itx2:jpr_itz2)%clgrid = 'V' ! and V-point
srcv(jpr_otx1:jpr_otz1)%laction = .TRUE. ! receive oce components on grid 1 only
srcv(jpr_itx1:jpr_itz2)%laction = .TRUE. ! receive ice components on grid 1 & 2
CASE default
CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_tau%clvgrd' )
END SELECT
IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' ) srcv(jpr_otx1:jpr_itz1)%nsgn = -1.
! ! Set grid and action
srcv(jpr_otx1:jpr_otz1)%laction = .TRUE. ! receive oce components on grid 1
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1
!
IF( TRIM( sn_rcv_tau%clvref ) == 'spherical' ) & ! spherical: 3rd component not received
& srcv( (/jpr_otz1, jpr_otz2, jpr_itz1, jpr_itz2/) )%laction = .FALSE.
!
IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' ) THEN ! already on local grid -> no need of the second grid
srcv(jpr_otx2:jpr_otz2)%laction = .FALSE.
srcv(jpr_itx2:jpr_itz2)%laction = .FALSE.
srcv(jpr_oty1)%clgrid = srcv(jpr_oty2)%clgrid ! not needed but cleaner...
srcv(jpr_ity1)%clgrid = srcv(jpr_ity2)%clgrid ! not needed but cleaner...
ENDIF
& srcv( (/jpr_otz1, jpr_itz1/) )%laction = .FALSE.
!
IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN ! 'oce and ice' case ocean stress on ocean mesh used
srcv(jpr_itx1:jpr_itz2)%laction = .FALSE. ! ice components not received
srcv(jpr_itx1)%clgrid = 'U' ! ocean stress used after its transformation
srcv(jpr_ity1)%clgrid = 'V' ! i.e. it is always at U- & V-points for i- & j-comp. resp.
IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN ! 'oce and ice' case ocean stress on ocean mesh used
srcv(jpr_itx1:jpr_itz1)%laction = .FALSE. ! ice components not received
ENDIF
ENDIF
......@@ -612,18 +549,18 @@ CONTAINS
ENDIF
srcv(jpr_topm:jpr_botm)%laction = .TRUE.
ENDIF
!! ! ! --------------------------- !
!! ! ! transmitted solar thru ice !
!! ! ! --------------------------- !
!! srcv(jpr_qtrice)%clname = 'OQtr'
!! IF( TRIM(sn_rcv_qtrice%cldes) == 'coupled' ) THEN
!! IF ( TRIM( sn_rcv_qtrice%clcat ) == 'yes' ) THEN
!! srcv(jpr_qtrice)%nct = nn_cats_cpl
!! ELSE
!! CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qtrice%clcat should always be set to yes currently' )
!! ENDIF
!! srcv(jpr_qtrice)%laction = .TRUE.
!! ENDIF
! ! --------------------------- !
! ! transmitted solar thru ice !
! ! --------------------------- !
srcv(jpr_qtrice)%clname = 'OQtr'
IF( TRIM(sn_rcv_qtrice%cldes) == 'coupled' ) THEN
IF ( TRIM( sn_rcv_qtrice%clcat ) == 'yes' ) THEN
srcv(jpr_qtrice)%nct = nn_cats_cpl
ELSE
CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qtrice%clcat should always be set to yes currently' )
ENDIF
srcv(jpr_qtrice)%laction = .TRUE.
ENDIF
! ! ------------------------- !
! ! ice skin temperature !
! ! ------------------------- !
......@@ -725,11 +662,10 @@ CONTAINS
srcv(:)%clgrid = 'T' ! force default definition in case of opa <-> sas coupling
srcv(:)%nsgn = 1. ! force default definition in case of opa <-> sas coupling
srcv( (/jpr_qsroce, jpr_qnsoce, jpr_oemp, jpr_sflx, jpr_fice, jpr_otx1, jpr_oty1, jpr_taum/) )%laction = .TRUE.
srcv(jpr_otx1)%clgrid = 'U' ! oce components given at U-point
srcv(jpr_oty1)%clgrid = 'V' ! and V-point
srcv(jpr_otx1)%clgrid = 'T' ! oce components given at T-point
srcv(jpr_oty1)%clgrid = 'T'
! Vectors: change of sign at north fold ONLY if on the local grid
srcv( (/jpr_otx1,jpr_oty1/) )%nsgn = -1.
sn_rcv_tau%clvgrd = 'U,V'
sn_rcv_tau%clvor = 'local grid'
sn_rcv_tau%clvref = 'spherical'
sn_rcv_emp%cldes = 'oce only'
......@@ -1162,7 +1098,7 @@ CONTAINS
!! ** Method : receive all fields from the atmosphere and transform
!! them into ocean surface boundary condition fields
!!
!! ** Action : update utau, vtau ocean stress at U,V grid
!! ** Action : update utau, vtau ocean stress at T-point
!! taum wind stress module at T-point
!! wndm wind speed module at T-point over free ocean or leads in presence of sea-ice
!! qns non solar heat fluxes including emp heat content (ocean only case)
......@@ -1221,39 +1157,20 @@ CONTAINS
IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN ! 2 components on the sphere
! ! (cartesian to spherical -> 3 to 2 components)
!
CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1), &
& srcv(jpr_otx1)%clgrid, ztx, zty )
CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1), 'T', ztx, zty )
frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:) ! overwrite 1st comp. on the 1st grid
frcv(jpr_oty1)%z3(:,:,1) = zty(:,:) ! overwrite 2nd comp. on the 1st grid
!
IF( srcv(jpr_otx2)%laction ) THEN
CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1), &
& srcv(jpr_otx2)%clgrid, ztx, zty )
frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:) ! overwrite 1st comp. on the 2nd grid
frcv(jpr_oty2)%z3(:,:,1) = zty(:,:) ! overwrite 2nd comp. on the 2nd grid
ENDIF
!
ENDIF
!
IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN ! 2 components oriented along the local grid
! ! (geographical to local grid -> rotate the components)
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )
IF( srcv(jpr_otx2)%laction ) THEN
CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )
ELSE
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )
ENDIF
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), 'T', 'en->i', ztx )
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), 'T', 'en->j', zty )
frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:) ! overwrite 1st component on the 1st grid
frcv(jpr_oty1)%z3(:,:,1) = zty(:,:) ! overwrite 2nd component on the 2nd grid
ENDIF
!
IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
DO_2D( 0, 0, 0, 0 ) ! T ==> (U,V)
frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
END_2D
CALL lbc_lnk( 'sbccpl', frcv(jpr_otx1)%z3(:,:,1), 'U', -1.0_wp, frcv(jpr_oty1)%z3(:,:,1), 'V', -1.0_wp )
ENDIF
llnewtx = .TRUE.
ELSE
llnewtx = .FALSE.
......@@ -1272,12 +1189,9 @@ CONTAINS
IF( .NOT. srcv(jpr_taum)%laction ) THEN ! compute wind stress module from its components if not received
! => need to be done only when otx1 was changed
IF( llnewtx ) THEN
DO_2D( 0, 0, 0, 0 )
zzx = frcv(jpr_otx1)%z3(ji-1,jj ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
zzy = frcv(jpr_oty1)%z3(ji ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
END_2D
CALL lbc_lnk( 'sbccpl', frcv(jpr_taum)%z3(:,:,1), 'T', 1.0_wp )
zzx = frcv(jpr_otx1)%z3(ji,jj,1)
zzy = frcv(jpr_oty1)%z3(ji,jj,1)
frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
llnewtau = .TRUE.
ELSE
llnewtau = .FALSE.
......@@ -1594,7 +1508,7 @@ CONTAINS
!! ** Action : return ptau_i, ptau_j, the stress over the ice
!!----------------------------------------------------------------------
REAL(wp), INTENT(inout), DIMENSION(:,:) :: p_taui ! i- & j-components of atmos-ice stress [N/m2]
REAL(wp), INTENT(inout), DIMENSION(:,:) :: p_tauj ! at I-point (B-grid) or U & V-point (C-grid)
REAL(wp), INTENT(inout), DIMENSION(:,:) :: p_tauj ! at T-point
!!
INTEGER :: ji, jj ! dummy loop indices
INTEGER :: itx ! index of taux over ice
......@@ -1616,28 +1530,16 @@ CONTAINS
!
IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN ! 2 components on the sphere
! ! (cartesian to spherical -> 3 to 2 components)
CALL geo2oce( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1), &
& srcv(jpr_itx1)%clgrid, ztx, zty )
CALL geo2oce( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1), 'T', ztx, zty )
frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:) ! overwrite 1st comp. on the 1st grid
frcv(jpr_ity1)%z3(:,:,1) = zty(:,:) ! overwrite 2nd comp. on the 1st grid
!
IF( srcv(jpr_itx2)%laction ) THEN
CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1), &
& srcv(jpr_itx2)%clgrid, ztx, zty )
frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:) ! overwrite 1st comp. on the 2nd grid
frcv(jpr_ity2)%z3(:,:,1) = zty(:,:) ! overwrite 2nd comp. on the 2nd grid
ENDIF
!
ENDIF
!
IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN ! 2 components oriented along the local grid
! ! (geographical to local grid -> rotate the components)
CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )
IF( srcv(jpr_itx2)%laction ) THEN
CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )
ELSE
CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty )
ENDIF
CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), 'T', 'en->i', ztx )
CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), 'T', 'en->j', zty )
frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:) ! overwrite 1st component on the 1st grid
frcv(jpr_ity1)%z3(:,:,1) = zty(:,:) ! overwrite 2nd component on the 1st grid
ENDIF
......@@ -1651,29 +1553,8 @@ CONTAINS
! ! ======================= !
! ! put on ice grid !
! ! ======================= !
!
! j+1 j -----V---F
! ice stress on ice velocity point ! |
! (C-grid ==>(U,V)) j | T U
! | |
! j j-1 -I-------|
! (for I) | |
! i-1 i i
! i i+1 (for I)
SELECT CASE ( srcv(jpr_itx1)%clgrid )
CASE( 'U' )
p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1) ! (U,V) ==> (U,V)
p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
CASE( 'T' )
DO_2D( 0, 0, 0, 0 ) ! T ==> (U,V)
! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology
zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) )
zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) )
p_taui(ji,jj) = zztmp1 * ( frcv(jpr_itx1)%z3(ji+1,jj ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
END_2D
CALL lbc_lnk( 'sbccpl', p_taui, 'U', -1., p_tauj, 'V', -1. )
END SELECT
p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)
p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
ENDIF
!
......@@ -1794,8 +1675,8 @@ CONTAINS
END SELECT
! --- evaporation over ice (kg/m2/s) --- !
IF (ln_scale_ice_flux) THEN ! typically met-office requirements
IF (sn_rcv_emp%clcat == 'yes') THEN
IF( ln_scale_ice_flux ) THEN ! typically met-office requirements
IF( sn_rcv_emp%clcat == 'yes' ) THEN
WHERE( a_i(:,:,:) > 1.e-10 ) ; zevap_ice(:,:,:) = frcv(jpr_ievp)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
ELSEWHERE ; zevap_ice(:,:,:) = 0._wp
END WHERE
......@@ -1812,7 +1693,7 @@ CONTAINS
ENDDO
ENDIF
ELSE
IF (sn_rcv_emp%clcat == 'yes') THEN
IF( sn_rcv_emp%clcat == 'yes' ) THEN
zevap_ice(:,:,1:jpl) = frcv(jpr_ievp)%z3(:,:,1:jpl)
WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
ELSEWHERE ; zevap_ice_total(:,:) = 0._wp
......@@ -1826,7 +1707,7 @@ CONTAINS
ENDIF
ENDIF
IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
IF( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
! For conservative case zemp_ice has not been defined yet. Do it now.
zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:)
ENDIF
......@@ -1926,13 +1807,13 @@ CONTAINS
& - zevap_ice_total(:,:) * picefr(:,:) ) * tmask(:,:,1) ) ! ice-free oce evap (cell average)
! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf
!! IF( srcv(jpr_rnf)%laction ) CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1) ) ! runoff
!! IF( srcv(jpr_isf)%laction ) CALL iom_put( 'iceshelf_cea', fwfisf(:,:) * tmask(:,:,1) ) ! iceshelf
!! IF( srcv(jpr_isf)%laction ) CALL iom_put( 'iceshelf_cea', fwfisf(:,:) * tmask(:,:,1) ) ! iceshelf
!
! ! ========================= !
SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) ) ! ice topmelt and botmelt !
! ! ========================= !
CASE ('coupled')
IF (ln_scale_ice_flux) THEN
CASE( 'coupled' )
IF( ln_scale_ice_flux ) THEN
WHERE( a_i(:,:,:) > 1.e-10_wp )
qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
......@@ -2213,28 +2094,29 @@ CONTAINS
!
ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN !== conduction flux as surface forcing ==!
!
!! SELECT CASE( TRIM( sn_rcv_qtrice%cldes ) )
!! !
!! ! ! ===> here we receive the qtr_ice_top array from the coupler
!! CASE ('coupled')
!! IF (ln_scale_ice_flux) THEN
!! WHERE( a_i(:,:,:) > 1.e-10_wp )
!! zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
!! ELSEWHERE
!! zqtr_ice_top(:,:,:) = 0.0_wp
!! ENDWHERE
!! ELSE
!! zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:)
!! ENDIF
!!
!! ! Add retrieved transmitted solar radiation onto the ice and total solar radiation
!! zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
!! zqsr_tot(:,:) = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
!!
!! ! if we are not getting this data from the coupler then assume zero (fully opaque ice)
!! CASE ('none')
zqtr_ice_top(:,:,:) = 0._wp
!! END SELECT
!
SELECT CASE( TRIM( sn_rcv_qtrice%cldes ) )
!
! ! ===> here we receive the qtr_ice_top array from the coupler
CASE ('coupled')
IF (ln_scale_ice_flux) THEN
WHERE( a_i(:,:,:) > 1.e-10_wp )
zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
ELSEWHERE
zqtr_ice_top(:,:,:) = 0.0_wp
ENDWHERE
ELSE
zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:)
ENDIF
! Add retrieved transmitted solar radiation onto the ice and total solar radiation
zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
zqsr_tot(:,:) = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
! if we are not getting this data from the coupler then assume zero (fully opaque ice)
CASE ('none')
zqtr_ice_top(:,:,:) = 0._wp
END SELECT
!
ENDIF
......@@ -2403,10 +2285,10 @@ CONTAINS
CASE( 'weighted ice' ) ;
SELECT CASE( sn_snd_alb%clcat )
CASE( 'yes' )
ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
CASE( 'no' )
WHERE( fr_i (:,:) > 0. )
ztmp1(:,:) = SUM ( alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 )
ztmp1(:,:) = SUM ( alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 )
ELSEWHERE
ztmp1(:,:) = 0.
END WHERE
......@@ -2572,17 +2454,18 @@ CONTAINS
IF( ssnd(jps_ocx1)%laction ) THEN ! Surface current !
! ! ------------------------- !
!
! j+1 j -----V---F
! surface velocity always sent from T point ! |
! j | T U
! | |
! j j-1 -I-------|
! (for I) | |
! i-1 i i
! i i+1 (for I)
! j -----V---F
! surface velocity always sent from T point ! |
! j | T U
! | |
! j-1 -I-------|
! | |
! i-1 i i
!!clem: make a new variable at T-point to replace uu and vv => uuT and vvT for instance
IF( nn_components == jp_iam_oce ) THEN
zotx1(:,:) = uu(:,:,1,Kmm)
zoty1(:,:) = vv(:,:,1,Kmm)
!!clem : should be demi sum, no? Or uuT and vvT
ELSE
SELECT CASE( TRIM( sn_snd_crt%cldes ) )
CASE( 'oce only' ) ! C-grid ==> T
......@@ -2652,42 +2535,42 @@ CONTAINS
! ! Surface current to waves !
! ! ------------------------- !
IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN
!
! j+1 j -----V---F
! surface velocity always sent from T point ! |
! j | T U
! | |
! j j-1 -I-------|
! (for I) | |
! i-1 i i
! i i+1 (for I)
SELECT CASE( TRIM( sn_snd_crtw%cldes ) )
CASE( 'oce only' ) ! C-grid ==> T
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj ,1,Kmm) )
zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji , jj-1,1,Kmm) )
END_2D
CASE( 'weighted oce and ice' ) ! Ocean and Ice on C-grid ==> T
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu (ji,jj,1,Kmm) + uu (ji-1,jj ,1,Kmm) ) * zfr_l(ji,jj)
zoty1(ji,jj) = 0.5 * ( vv (ji,jj,1,Kmm) + vv (ji ,jj-1,1,Kmm) ) * zfr_l(ji,jj)
zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * fr_i(ji,jj)
zity1(ji,jj) = 0.5 * ( v_ice(ji,jj ) + v_ice(ji ,jj-1 ) ) * fr_i(ji,jj)
END_2D
CALL lbc_lnk( 'sbccpl', zitx1, 'T', -1.0_wp, zity1, 'T', -1.0_wp )
CASE( 'mixed oce-ice' ) ! Ocean and Ice on C-grid ==> T
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu (ji,jj,1,Kmm) + uu (ji-1,jj ,1,Kmm) ) * zfr_l(ji,jj) &
& + 0.5 * ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * fr_i(ji,jj)
zoty1(ji,jj) = 0.5 * ( vv (ji,jj,1,Kmm) + vv (ji ,jj-1,1,Kmm) ) * zfr_l(ji,jj) &
& + 0.5 * ( v_ice(ji,jj ) + v_ice(ji ,jj-1 ) ) * fr_i(ji,jj)
END_2D
END SELECT
!
! j -----V---F
! surface velocity always sent from T point ! |
! j | T U
! | |
! j-1 -I-------|
! | |
! i-1 i i
!!clem: make a new variable at T-point to replace uu and vv => uuT and vvT for instance
SELECT CASE( TRIM( sn_snd_crtw%cldes ) )
CASE( 'oce only' ) ! C-grid ==> T
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj ,1,Kmm) )
zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji , jj-1,1,Kmm) )
END_2D
CASE( 'weighted oce and ice' ) ! Ocean and Ice on C-grid ==> T
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu (ji,jj,1,Kmm) + uu (ji-1,jj ,1,Kmm) ) * zfr_l(ji,jj)
zoty1(ji,jj) = 0.5 * ( vv (ji,jj,1,Kmm) + vv (ji ,jj-1,1,Kmm) ) * zfr_l(ji,jj)
zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * fr_i(ji,jj)
zity1(ji,jj) = 0.5 * ( v_ice(ji,jj ) + v_ice(ji ,jj-1 ) ) * fr_i(ji,jj)
END_2D
CALL lbc_lnk( 'sbccpl', zitx1, 'T', -1.0_wp, zity1, 'T', -1.0_wp )
CASE( 'mixed oce-ice' ) ! Ocean and Ice on C-grid ==> T
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu (ji,jj,1,Kmm) + uu (ji-1,jj ,1,Kmm) ) * zfr_l(ji,jj) &
& + 0.5 * ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * fr_i(ji,jj)
zoty1(ji,jj) = 0.5 * ( vv (ji,jj,1,Kmm) + vv (ji ,jj-1,1,Kmm) ) * zfr_l(ji,jj) &
& + 0.5 * ( v_ice(ji,jj ) + v_ice(ji ,jj-1 ) ) * fr_i(ji,jj)
END_2D
END SELECT
CALL lbc_lnk( 'sbccpl', zotx1, ssnd(jps_ocxw)%clgrid, -1.0_wp, zoty1, ssnd(jps_ocyw)%clgrid, -1.0_wp )
!
!
IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN ! Rotation of the components
! ! Ocean component
! ! Ocean component
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 ) ! 1st component
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 ) ! 2nd component
zotx1(:,:) = ztmp1(:,:) ! overwrite the components
......@@ -2700,18 +2583,18 @@ CONTAINS
ENDIF
ENDIF
!
! ! spherical coordinates to cartesian -> 2 components to 3 components
! IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN
! ztmp1(:,:) = zotx1(:,:) ! ocean currents
! ztmp2(:,:) = zoty1(:,:)
! CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
! !
! IF( ssnd(jps_ivx1)%laction ) THEN ! ice velocities
! ztmp1(:,:) = zitx1(:,:)
! ztmp1(:,:) = zity1(:,:)
! CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
! ENDIF
! ENDIF
! ! spherical coordinates to cartesian -> 2 components to 3 components
! IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN
! ztmp1(:,:) = zotx1(:,:) ! ocean currents
! ztmp2(:,:) = zoty1(:,:)
! CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
! !
! IF( ssnd(jps_ivx1)%laction ) THEN ! ice velocities
! ztmp1(:,:) = zitx1(:,:)
! ztmp1(:,:) = zity1(:,:)
! CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
! ENDIF
! ENDIF
!
IF( ssnd(jps_ocxw)%laction ) CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info ) ! ocean x current 1st grid
IF( ssnd(jps_ocyw)%laction ) CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info ) ! ocean y current 1st grid
......
src/OCE/SBC/sbcflx.F90
View file @ ee6ce8c5
......@@ -119,8 +119,8 @@ CONTAINS
END DO
! ! fill sf with slf_i and control print
CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
sf(jp_utau)%cltype = 'U' ; sf(jp_utau)%zsgn = -1._wp ! vector field at U point: overwrite default definition of cltype and zsgn
sf(jp_vtau)%cltype = 'V' ; sf(jp_vtau)%zsgn = -1._wp ! vector field at V point: overwrite default definition of cltype and zsgn
sf(jp_utau)%cltype = 'T' ; sf(jp_utau)%zsgn = -1._wp ! vector field at T point: overwrite default definition of cltype and zsgn
sf(jp_vtau)%cltype = 'T' ; sf(jp_vtau)%zsgn = -1._wp ! vector field at T point: overwrite default definition of cltype and zsgn
!
ENDIF
......@@ -147,8 +147,8 @@ CONTAINS
ENDIF
#endif
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! set the ocean fluxes from read fields
utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) * umask(ji,jj,1)
vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) * vmask(ji,jj,1)
utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) * tmask(ji,jj,1)
vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) * tmask(ji,jj,1)
qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) * tmask(ji,jj,1)
!!sfx (ji,jj) = sf(jp_sfx )%fnow(ji,jj,1) * tmask(ji,jj,1)
......@@ -170,19 +170,15 @@ CONTAINS
ENDIF
!
ENDIF
! ! module of wind stress and wind speed at T-point
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
!
! module of wind stress and wind speed at T-point
zcoef = 1. / ( zrhoa * zcdrag )
DO_2D( 0, 0, 0, 0 )
ztx = ( utau(ji-1,jj ) + utau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( umask(ji-1,jj ,1), umask(ji,jj,1) ) )
zty = ( vtau(ji ,jj-1) + vtau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( vmask(ji ,jj-1,1), vmask(ji,jj,1) ) )
zmod = SQRT( ztx * ztx + zty * zty ) * tmask(ji,jj,1)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zmod = SQRT( utau(ji,jj) * utau(ji,jj) + vtau(ji,jj) * vtau(ji,jj) ) * tmask(ji,jj,1)
taum(ji,jj) = zmod
wndm(ji,jj) = SQRT( zmod * zcoef ) !!clem: not used?
END_2D
!
CALL lbc_lnk( 'sbcflx', taum, 'T', 1._wp, wndm, 'T', 1._wp )
!
END SUBROUTINE sbc_flx
!!======================================================================
......
src/OCE/SBC/sbcmod.F90
View file @ ee6ce8c5
......@@ -158,8 +158,8 @@ CONTAINS
!
IF( .NOT.ln_usr ) THEN ! the model calendar needs some specificities (except in user defined case)
IF( MOD( rday , rn_Dt ) /= 0. ) CALL ctl_stop( 'the time step must devide the number of second of in a day' )
IF( MOD( rday , 2. ) /= 0. ) CALL ctl_stop( 'the number of second of in a day must be an even number' )
IF( MOD( rn_Dt , 2. ) /= 0. ) CALL ctl_stop( 'the time step (in second) must be an even number' )
IF( MOD( rday , 2. ) /= 0. ) CALL ctl_stop( 'the number of second of in a day must be an even number' )
IF( MOD( rn_Dt, 2. ) /= 0. ) CALL ctl_stop( 'the time step (in second) must be an even number' )
ENDIF
! !** check option consistency
!
......@@ -395,16 +395,16 @@ CONTAINS
! ! ---------------------------------------- !
IF( kt /= nit000 ) THEN ! Swap of forcing fields !
! ! ---------------------------------------- !
utau_b(:,:) = utau(:,:) ! Swap the ocean forcing fields
vtau_b(:,:) = vtau(:,:) ! (except at nit000 where before fields
qns_b (:,:) = qns (:,:) ! are set at the end of the routine)
emp_b (:,:) = emp (:,:)
sfx_b (:,:) = sfx (:,:)
utau_b(:,:) = utauU(:,:) ! Swap the ocean forcing fields
vtau_b(:,:) = vtauV(:,:) ! (except at nit000 where before fields
qns_b (:,:) = qns (:,:) ! are set at the end of the routine)
emp_b (:,:) = emp (:,:)
sfx_b (:,:) = sfx (:,:)
IF( ln_rnf ) THEN
rnf_b (:,: ) = rnf (:,: )
rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
ENDIF
!
!
ENDIF
! ! ---------------------------------------- !
! ! forcing field computation !
......@@ -420,60 +420,62 @@ CONTAINS
!
SELECT CASE( nsbc ) ! Compute ocean surface boundary condition
! ! (i.e. utau,vtau, qns, qsr, emp, sfx)
CASE( jp_usr ) ; CALL usrdef_sbc_oce( kt, Kbb ) ! user defined formulation
CASE( jp_flx ) ; CALL sbc_flx ( kt ) ! flux formulation
CASE( jp_usr ) ; CALL usrdef_sbc_oce( kt, Kbb ) ! user defined formulation
CASE( jp_flx ) ; CALL sbc_flx ( kt ) ! flux formulation
CASE( jp_blk )
IF( ll_sas ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-SAS coupling: SAS receiving fields from OCE
!!!!!!!!!!! ATTENTION:ln_wave is not only used for oasis coupling !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( ln_wave ) THEN
IF ( lk_oasis ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-wave coupling
IF ( lk_oasis ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-wave coupling
CALL sbc_wave ( kt, Kmm )
ENDIF
CALL sbc_blk ( kt ) ! bulk formulation for the ocean
CALL sbc_blk ( kt ) ! bulk formulation for the ocean
!
CASE( jp_abl )
IF( ll_sas ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-SAS coupling: SAS receiving fields from OCE
CALL sbc_abl ( kt ) ! ABL formulation for the ocean
CALL sbc_abl ( kt ) ! ABL formulation for the ocean
!
CASE( jp_purecpl ) ; CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! pure coupled formulation
CASE( jp_none )
IF( ll_opa ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-SAS coupling: OCE receiving fields from SAS
IF( ll_opa ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-SAS coupling: OCE receiving fields from SAS
END SELECT
!
IF( ln_mixcpl ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! forced-coupled mixed formulation after forcing
IF( ln_mixcpl ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! forced-coupled mixed formulation after forcing
!
IF( ln_wave .AND. ln_tauoc ) THEN ! Wave stress reduction
DO_2D( 0, 0, 0, 0)
utau(ji,jj) = utau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji+1,jj) ) * 0.5_wp
vtau(ji,jj) = vtau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji,jj+1) ) * 0.5_wp
END_2D
!
CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
!
taum(:,:) = taum(:,:)*tauoc_wave(:,:)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
utau(ji,jj) = utau(ji,jj) * tauoc_wave(ji,jj)
vtau(ji,jj) = vtau(ji,jj) * tauoc_wave(ji,jj)
taum(ji,jj) = taum(ji,jj) * tauoc_wave(ji,jj)
END_2D
!
IF( kt == nit000 ) CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.', &
& 'If not requested select ln_tauoc=.false.' )
!
ELSEIF( ln_wave .AND. ln_taw ) THEN ! Wave stress reduction
utau(:,:) = utau(:,:) - tawx(:,:) + twox(:,:)
vtau(:,:) = vtau(:,:) - tawy(:,:) + twoy(:,:)
CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
!
DO_2D( 0, 0, 0, 0)
taum(ji,jj) = sqrt((.5*(utau(ji-1,jj)+utau(ji,jj)))**2 + (.5*(vtau(ji,jj-1)+vtau(ji,jj)))**2)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
utau(ji,jj) = utau(ji,jj) - tawx(ji,jj) + twox(ji,jj)
vtau(ji,jj) = vtau(ji,jj) - tawy(ji,jj) + twoy(ji,jj)
taum(ji,jj) = SQRT( utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj) )
END_2D
!
IF( kt == nit000 ) CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.', &
& 'If not requested select ln_taw=.false.' )
!
ENDIF
CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1. )
!
!
IF( ln_icebergs ) THEN ! save pure stresses (with no ice-ocean stress) for use by icebergs
utau_icb(:,:) = utau(:,:) ; vtau_icb(:,:) = vtau(:,:)
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
! and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
DO_2D( 0, 0, 0, 0 )
utau_icb(ji,jj) = 0.5_wp * ( utau(ji,jj) + utau(ji+1,jj) ) * &
& ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
vtau_icb(ji,jj) = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj+1) ) * &
& ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
END_2D
CALL lbc_lnk( 'sbcmod', utau_icb, 'U', -1.0_wp, vtau_icb, 'V', -1.0_wp )
ENDIF
!
! !== Misc. Options ==!
......@@ -507,9 +509,6 @@ CONTAINS
! Should not be run if ln_diurnal_only
IF( l_sbc_clo ) CALL sbc_clo( kt )
!!$!RBbug do not understand why see ticket 667
!!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
!!$ CALL lbc_lnk( 'sbcmod', emp, 'T', 1.0_wp )
IF( ll_wd ) THEN ! If near WAD point limit the flux for now
zthscl = atanh(rn_wd_sbcfra) ! taper frac default is .999
zwdht(:,:) = ssh(:,:,Kmm) + ht_0(:,:) - rn_wdmin1 ! do this calc of water
......@@ -534,6 +533,16 @@ CONTAINS
emp (:,:) = emp(:,:) * zwght(:,:)
END WHERE
ENDIF
! --- calculate utau and vtau on U,V-points --- !
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
! and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
utauU (ji,jj) = 0.5_wp * ( utau(ji,jj) + utau(ji+1,jj) ) * &
& ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
vtauV (ji,jj) = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj+1) ) * &
& ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
END_2D
IF( nn_hls == 1 ) CALL lbc_lnk( 'sbcmod', utauU, 'U', -1.0_wp, vtauV, 'V', -1.0_wp )
!
IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 !
! ! ---------------------------------------- !
......@@ -543,27 +552,28 @@ CONTAINS
IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN !* MLF: Restart: read in restart file
#endif
IF(lwp) WRITE(numout,*) ' nit000-1 surface forcing fields read in the restart file'
CALL iom_get( numror, jpdom_auto, 'utau_b', utau_b ) ! i-stress
CALL iom_get( numror, jpdom_auto, 'vtau_b', vtau_b ) ! j-stress
CALL iom_get( numror, jpdom_auto, 'qns_b', qns_b ) ! non solar heat flux
CALL iom_get( numror, jpdom_auto, 'emp_b', emp_b ) ! freshwater flux
CALL iom_get( numror, jpdom_auto, 'utau_b', utau_b, cd_type = 'U', psgn = -1._wp ) ! i-stress
CALL iom_get( numror, jpdom_auto, 'vtau_b', vtau_b, cd_type = 'V', psgn = -1._wp ) ! j-stress
CALL iom_get( numror, jpdom_auto, 'qns_b', qns_b, cd_type = 'T', psgn = 1._wp ) ! non solar heat flux
CALL iom_get( numror, jpdom_auto, 'emp_b', emp_b, cd_type = 'T', psgn = 1._wp ) ! freshwater flux
! NB: The 3D heat content due to qsr forcing (qsr_hc_b) is treated in traqsr
! To ensure restart capability with 3.3x/3.4 restart files !! to be removed in v3.6
IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN
CALL iom_get( numror, jpdom_auto, 'sfx_b', sfx_b ) ! before salt flux (T-point)
CALL iom_get( numror, jpdom_auto, 'sfx_b', sfx_b, cd_type = 'T', psgn = 1._wp ) ! before salt flux (T-point)
ELSE
sfx_b (:,:) = sfx(:,:)
ENDIF
ELSE !* no restart: set from nit000 values
IF(lwp) WRITE(numout,*) ' nit000-1 surface forcing fields set to nit000'
utau_b(:,:) = utau(:,:)
vtau_b(:,:) = vtau(:,:)
utau_b(:,:) = utauU(:,:)
vtau_b(:,:) = vtauV(:,:)
qns_b (:,:) = qns (:,:)
emp_b (:,:) = emp (:,:)
sfx_b (:,:) = sfx (:,:)
ENDIF
ENDIF
!
!
#if defined key_RK3
! ! ---------------------------------------- !
IF( lrst_oce .AND. lk_SWE ) THEN ! RK3: Write in the ocean restart file !
......@@ -578,13 +588,13 @@ CONTAINS
IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ', &
& 'at it= ', kt,' date= ', ndastp
IF(lwp) WRITE(numout,*) '~~~~'
CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )
CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau )
CALL iom_rstput( kt, nitrst, numrow, 'qns_b' , qns )
CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utauU )
CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtauV )
CALL iom_rstput( kt, nitrst, numrow, 'qns_b' , qns )
! The 3D heat content due to qsr forcing is treated in traqsr
! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b' , qsr )
CALL iom_rstput( kt, nitrst, numrow, 'emp_b' , emp )
CALL iom_rstput( kt, nitrst, numrow, 'sfx_b' , sfx )
CALL iom_rstput( kt, nitrst, numrow, 'emp_b' , emp )
CALL iom_rstput( kt, nitrst, numrow, 'sfx_b' , sfx )
ENDIF
! ! ---------------------------------------- !
! ! Outputs and control print !
......@@ -628,8 +638,8 @@ CONTAINS
CALL prt_ctl(tab3d_1=tmask , clinfo1=' tmask - : ', mask1=tmask, kdim=jpk )
CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' sst - : ', mask1=tmask, kdim=1 )
CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sss - : ', mask1=tmask, kdim=1 )
CALL prt_ctl(tab2d_1=utau , clinfo1=' utau - : ', mask1=umask, &
& tab2d_2=vtau , clinfo2=' vtau - : ', mask2=vmask )
CALL prt_ctl(tab2d_1=utau , clinfo1=' utau - : ', mask1=tmask, &
& tab2d_2=vtau , clinfo2=' vtau - : ', mask2=tmask )
ENDIF
IF( kt == nitend ) CALL sbc_final ! Close down surface module if necessary
......
src/OCE/TRA/eosbn2.F90
View file @ ee6ce8c5
......@@ -1336,10 +1336,10 @@ CONTAINS
!! pab_pe(:,:,:,jp_tem) is alpha_pe
!! pab_pe(:,:,:,jp_sal) is beta_pe
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab_pe ! alpha_pe and beta_pe
REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: ppen ! potential energy anomaly
INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity
REAL(wp), DIMENSION(:,:,:,:), INTENT( out) :: pab_pe ! alpha_pe and beta_pe
REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: ppen ! potential energy anomaly
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs , ztm ! local scalars
......@@ -1352,7 +1352,7 @@ CONTAINS
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
!
zh = gdept(ji,jj,jk,Kmm) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
......@@ -1411,9 +1411,9 @@ CONTAINS
!
CASE( np_seos ) !== Vallis (2006) simplified EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zt = pts(ji,jj,jk,jp_tem) - 10._wp ! temperature anomaly (t-T0)
zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point
ztm = tmask(ji,jj,jk) ! tmask
zn = 0.5_wp * zh * r1_rho0 * ztm
......
src/OCE/TRA/traadv.F90
View file @ ee6ce8c5
......@@ -10,6 +10,7 @@ MODULE traadv
!! - ! 2014-12 (G. Madec) suppression of cross land advection option
!! 3.6 ! 2015-06 (E. Clementi) Addition of Stokes drift in case of wave coupling
!! 4.5 ! 2021-04 (G. Madec, S. Techene) add advective velocities as optional arguments
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -191,11 +192,19 @@ CONTAINS
SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==!
!
CASE ( np_CEN ) ! Centered scheme : 2nd / 4th order
CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
IF( nn_hls == 1 ) THEN
CALL tra_adv_cen_hls1( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
ELSE
CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
ENDIF
CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order
CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )
CASE ( np_MUS ) ! MUSCL
IF( nn_hls == 1 ) THEN
CALL tra_adv_mus_hls1( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
ELSE
CALL tra_adv_mus( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
END IF
CASE ( np_UBS ) ! UBS
CALL tra_adv_ubs ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v )
CASE ( np_QCK ) ! QUICKEST
......
src/OCE/TRA/traadv_cen.F90
View file @ ee6ce8c5
......@@ -4,6 +4,7 @@ MODULE traadv_cen
!! Ocean tracers: advective trend (2nd/4th order centered)
!!======================================================================
!! History : 3.7 ! 2014-05 (G. Madec) original code
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -29,7 +30,8 @@ MODULE traadv_cen
IMPLICIT NONE
PRIVATE
PUBLIC tra_adv_cen ! called by traadv.F90
PUBLIC tra_adv_cen ! called by traadv.F90
PUBLIC tra_adv_cen_hls1 ! called by traadv.F90
REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6
......@@ -47,7 +49,321 @@ MODULE traadv_cen
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE tra_adv_cen_test( kt, kit000, cdtype, pU, pV, pW, &
& Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_cen ***
!!
!! ** Purpose : Compute the now trend due to the advection of tracers
!! and add it to the general trend of passive tracer equations.
!!
!! ** Method : The advection is evaluated by a 2nd or 4th order scheme
!! using now fields (leap-frog scheme).
!! kn_cen_h = 2 ==>> 2nd order centered scheme on the horizontal
!! = 4 ==>> 4th order - - - -
!! kn_cen_v = 2 ==>> 2nd order centered scheme on the vertical
!! = 4 ==>> 4th order COMPACT scheme - -
!!
!! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
!! - poleward advective heat and salt transport (l_diaptr=T)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
INTEGER , INTENT(in ) :: kit000 ! first time step index
CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
INTEGER , INTENT(in ) :: kjpt ! number of tracers
INTEGER , INTENT(in ) :: kn_cen_h ! =2/4 (2nd or 4th order scheme)
INTEGER , INTENT(in ) :: kn_cen_v ! =2/4 (2nd or 4th order scheme)
! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
INTEGER :: ierr ! local integer
REAL(wp) :: zC2t_u, zC4t_u ! local scalars
REAL(wp) :: zC2t_v, zC4t_v ! - -
REAL(wp) :: zftw_kp1
REAL(wp), DIMENSION(A2D(1)) :: zft_u, zft_v !, zft_w
REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zdt_u, zdt_v
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: ztw
!!----------------------------------------------------------------------
!
#if defined key_loop_fusion
CALL tra_adv_cen_lf ( kt, nit000, cdtype, pU, pV, pW, Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
#else
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == kit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
ENDIF
! ! set local switches
l_trd = .FALSE.
l_hst = .FALSE.
l_ptr = .FALSE.
IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
ENDIF
!
! !* 2nd order centered
DO jn = 1, kjpt !== loop over the tracers ==!
!
DO jk = 1, jpkm1
!
DO_2D( 1, 0, 1, 0 ) ! Horizontal fluxes at layer jk
zft_u(ji,jj) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) )
zft_v(ji,jj) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) )
END_2D
!
DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
& + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
END DO
!
#define zft_w zft_u
IF( ln_linssh ) THEN !* top value (linear free surf. only as zwz is multiplied by wmask)
DO_2D( 0, 0, 0, 0 )
zft_w(ji,jj) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
zft_w(ji,jj) = 0._wp
END_2D
ENDIF
DO jk = 1, jpk-2
DO_2D( 0, 0, 0, 0 ) ! Vertical fluxes
zftw_kp1 = 0.5 * pW(ji,jj,jk+1) * ( pt(ji,jj,jk+1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk+1)
!
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_w(ji,jj) - zftw_kp1 ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
zft_w(ji,jj) = zftw_kp1
END_2D
END DO
jk = jpkm1 ! bottom vertical flux set to zero for all tracers
DO_2D( 0, 0, 0, 0 )
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zft_w(ji,jj) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
!
END DO
!
!
#undef zft_w
! ! trend diagnostics
!!gm + !!st to be done with the whole rewritting of trd
!! trd routine arguments MUST be changed adding jk and zwx, zwy in 2D
!! IF( l_trd ) THEN
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
!! ENDIF
!! ! ! "Poleward" heat and salt transports
!! IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
!! ! ! heat and salt transport
!! IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
!
!
!
#endif
END SUBROUTINE tra_adv_cen_test
SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pU, pV, pW, &
& Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_cen ***
!!
!! ** Purpose : Compute the now trend due to the advection of tracers
!! and add it to the general trend of passive tracer equations.
!!
!! ** Method : The advection is evaluated by a 2nd or 4th order scheme
!! using now fields (leap-frog scheme).
!! kn_cen_h = 2 ==>> 2nd order centered scheme on the horizontal
!! = 4 ==>> 4th order - - - -
!! kn_cen_v = 2 ==>> 2nd order centered scheme on the vertical
!! = 4 ==>> 4th order COMPACT scheme - -
!!
!! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
!! - poleward advective heat and salt transport (l_diaptr=T)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
INTEGER , INTENT(in ) :: kit000 ! first time step index
CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
INTEGER , INTENT(in ) :: kjpt ! number of tracers
INTEGER , INTENT(in ) :: kn_cen_h ! =2/4 (2nd or 4th order scheme)
INTEGER , INTENT(in ) :: kn_cen_v ! =2/4 (2nd or 4th order scheme)
! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
INTEGER :: ierr ! local integer
REAL(wp) :: zC2t_u, zC4t_u ! local scalars
REAL(wp) :: zC2t_v, zC4t_v ! - -
REAL(wp) :: zftw_kp1
REAL(wp), DIMENSION(A2D(1)) :: zft_u, zft_v
REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zdt_u, zdt_v
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: ztw
!!----------------------------------------------------------------------
!
#if defined key_loop_fusion
CALL tra_adv_cen_lf ( kt, nit000, cdtype, pU, pV, pW, Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
#else
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == kit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
ENDIF
! ! set local switches
l_trd = .FALSE.
l_hst = .FALSE.
l_ptr = .FALSE.
IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
ENDIF
!
IF( kn_cen_h == 4 ) ALLOCATE( zdt_u(A2D(2)) , zdt_v(A2D(2)) ) ! horizontal 4th order only
IF( kn_cen_v == 4 ) ALLOCATE( ztw(A2D(nn_hls),jpk) ) ! vertical 4th order only
!
DO jn = 1, kjpt !== loop over the tracers ==!
!
SELECT CASE( kn_cen_h ) !-- Horizontal divergence of advective fluxes --!
!
!!st limitation : does not take into acccount iceshelf specificity
!! in case of linssh
CASE( 2 ) !* 2nd order centered
DO jk = 1, jpkm1
!
DO_2D( 1, 0, 1, 0 ) ! Horizontal fluxes at layer jk
zft_u(ji,jj) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) )
zft_v(ji,jj) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) )
END_2D
!
DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
& + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
END DO
!
CASE( 4 ) !* 4th order centered
DO jk = 1, jpkm1
DO_2D( 2, 1, 2, 1 ) ! masked gradient
zdt_u(ji,jj) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
zdt_v(ji,jj) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
END_2D
!
DO_2D( 1, 0, 1, 0 ) ! Horizontal advective fluxes
zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! C2 interpolation of T at u- & v-points (x2)
zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm)
! ! C4 interpolation of T at u- & v-points (x2)
zC4t_u = zC2t_u + r1_6 * ( zdt_u(ji-1,jj ) - zdt_u(ji+1,jj ) )
zC4t_v = zC2t_v + r1_6 * ( zdt_v(ji ,jj-1) - zdt_v(ji ,jj+1) )
! ! C4 fluxes
zft_u(ji,jj) = 0.5_wp * pU(ji,jj,jk) * zC4t_u
zft_v(ji,jj) = 0.5_wp * pV(ji,jj,jk) * zC4t_v
END_2D
!
DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
& + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
END DO
!
CASE DEFAULT
CALL ctl_stop( 'traadv_cen: wrong value for nn_cen' )
END SELECT
!
#define zft_w zft_u
!
IF( ln_linssh ) THEN !* top value (linear free surf. only as zwz is multiplied by wmask)
DO_2D( 0, 0, 0, 0 )
zft_w(ji,jj) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
zft_w(ji,jj) = 0._wp
END_2D
ENDIF
!
SELECT CASE( kn_cen_v ) !-- Vertical divergence of advective fluxes --! (interior)
!
CASE( 2 ) !* 2nd order centered
DO jk = 1, jpk-2
DO_2D( 0, 0, 0, 0 ) ! Vertical fluxes
zftw_kp1 = 0.5 * pW(ji,jj,jk+1) * ( pt(ji,jj,jk+1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk+1)
!
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_w(ji,jj) - zftw_kp1 ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
zft_w(ji,jj) = zftw_kp1
END_2D
END DO
jk = jpkm1 ! bottom vertical flux set to zero for all tracers
DO_2D( 0, 0, 0, 0 )
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zft_w(ji,jj) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
!
CASE( 4 ) !* 4th order compact
CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! ztw = interpolated value of T at w-point
!
DO jk = 1, jpk-2
!
DO_2D( 0, 0, 0, 0 )
zftw_kp1 = pW(ji,jj,jk+1) * ztw(ji,jj,jk+1) * wmask(ji,jj,jk+1)
! ! Divergence of advective fluxes
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_w(ji,jj) - zftw_kp1 ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
! ! update
zft_w(ji,jj) = zftw_kp1
END_2D
!
END DO
!
jk = jpkm1 ! bottom vertical flux set to zero for all tracers
DO_2D( 0, 0, 0, 0 )
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zft_w(ji,jj) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
!
END SELECT
!
#undef zft_w
! ! trend diagnostics
!!gm + !!st to be done with the whole rewritting of trd
!! trd routine arguments MUST be changed adding jk and zwx, zwy in 2D
!! IF( l_trd ) THEN
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
!! ENDIF
!! ! ! "Poleward" heat and salt transports
!! IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
!! ! ! heat and salt transport
!! IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
!
END DO
!
IF( kn_cen_h == 4 ) DEALLOCATE( zdt_u , zdt_v ) ! horizontal 4th order only
IF( kn_cen_v == 4 ) DEALLOCATE( ztw ) ! vertical 4th order only
!
#endif
END SUBROUTINE tra_adv_cen
SUBROUTINE tra_adv_cen_hls1( kt, kit000, cdtype, pU, pV, pW, &
& Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_cen ***
......@@ -141,6 +457,12 @@ CONTAINS
CASE DEFAULT
CALL ctl_stop( 'traadv_cen: wrong value for nn_cen' )
END SELECT
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --!
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
& - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
!
SELECT CASE( kn_cen_v ) !-- Vertical fluxes --! (interior)
!
......@@ -171,11 +493,17 @@ CONTAINS
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --!
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
& - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
& - ( zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
!
!!st DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --!
!!st pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
!!st & - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
!!st & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
!!st & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
!!st & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
!!st END_3D
! ! trend diagnostics
IF( l_trd ) THEN
CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
......@@ -190,7 +518,7 @@ CONTAINS
END DO
!
#endif
END SUBROUTINE tra_adv_cen
END SUBROUTINE tra_adv_cen_hls1
!!======================================================================
END MODULE traadv_cen
src/OCE/TRA/traadv_mus.F90
View file @ ee6ce8c5
......@@ -9,6 +9,7 @@ MODULE traadv_mus
!! 3.2 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
!! 3.4 ! 2012-06 (P. Oddo, M. Vichi) include the upstream where needed
!! 3.7 ! 2015-09 (G. Madec) add the ice-shelf cavities boundary condition
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -34,7 +35,9 @@ MODULE traadv_mus
IMPLICIT NONE
PRIVATE
PUBLIC tra_adv_mus ! routine called by traadv.F90
PUBLIC tra_adv_mus ! routine called by traadv.F90
PUBLIC tra_adv_mus_hls1 ! routine called by traadv.F90
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
! ! and in closed seas (orca 2 and 1 configurations)
......@@ -55,6 +58,217 @@ MODULE traadv_mus
CONTAINS
SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pU, pV, pW, &
& Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_mus ***
!!
!! ** Purpose : Compute the now trend due to total advection of tracers
!! using a MUSCL scheme (Monotone Upstream-centered Scheme for
!! Conservation Laws) and add it to the general tracer trend.
!!
!! ** Method : MUSCL scheme plus centered scheme at ocean boundaries
!! ld_msc_ups=T :
!!
!! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
!! - poleward advective heat and salt transport (ln_diaptr=T)
!!
!! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
!! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices
INTEGER , INTENT(in ) :: kit000 ! first time step index
CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
INTEGER , INTENT(in ) :: kjpt ! number of tracers
LOGICAL , INTENT(in ) :: ld_msc_ups ! use upstream scheme within muscl
REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
INTEGER :: ierr ! local integer
INTEGER :: ik
REAL(wp) :: zu, z0u, zzslpx, zzwx, zw , zalpha ! local scalars
REAL(wp) :: zv, z0v, zzslpy, zzwy, z0w ! - -
REAL(wp) :: zdzt_kp2, zslpz_kp1, zfW_kp1
REAL(wp), DIMENSION(A2D(2)) :: zdxt, zslpx, zwx ! 2D workspace
REAL(wp), DIMENSION(A2D(2)) :: zdyt, zslpy, zwy ! - -
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == kit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
IF(lwp) WRITE(numout,*) ' : mixed up-stream ', ld_msc_ups
IF(lwp) WRITE(numout,*) '~~~~~~~'
IF(lwp) WRITE(numout,*)
!
! Upstream / MUSCL scheme indicator
!
ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
xind(:,:,:) = 1._wp ! set equal to 1 where up-stream is not needed
IF( ld_msc_ups ) THEN ! define the upstream indicator (if asked)
DO jk = 1, jpkm1
xind(:,:,jk) = 1._wp & ! =>1 where up-stream is not needed
& - rnfmsk(:,:) * rnfmsk_z(jk) * tmask(:,:,jk) ! =>0 near runoff mouths (& closed sea outflows)
END DO
ENDIF
!
ENDIF
!
l_trd = .FALSE.
l_hst = .FALSE.
l_ptr = .FALSE.
IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
IF( l_diaptr .AND. cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
IF( l_iom .AND. cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
ENDIF
!
!
DO jn = 1, kjpt !== loop over the tracers ==!
!
DO jk = 1, jpkm1
! !* Horizontal advective fluxes
!
! !-- first guess of the slopes
DO_2D( 2, 1, 2, 1 )
zdxt(ji,jj) = ( pt(ji+1,jj ,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) * umask(ji,jj,jk)
zdyt(ji,jj) = ( pt(ji ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) * vmask(ji,jj,jk)
END_2D
! !-- Slopes of tracer
DO_2D( 1, 1, 1, 1 )
! ! 1/2 Slopes at T-point (set to 0 if adjectent slopes are of opposite sign)
zzslpx = ( zdxt(ji,jj) + zdxt(ji-1,jj ) ) &
& * ( 0.25_wp + SIGN( 0.25_wp, zdxt(ji,jj) * zdxt(ji-1,jj ) ) )
zzslpy = ( zdyt(ji,jj) + zdyt(ji ,jj-1) ) &
& * ( 0.25_wp + SIGN( 0.25_wp, zdyt(ji,jj) * zdyt(ji ,jj-1) ) )
! ! Slopes limitation
zslpx(ji,jj) = SIGN( 1.0_wp, zzslpx ) * MIN( ABS( zzslpx ), &
& 2._wp*ABS( zdxt (ji-1,jj) ), &
& 2._wp*ABS( zdxt (ji ,jj) ) )
zslpy(ji,jj) = SIGN( 1.0_wp, zzslpy ) * MIN( ABS( zzslpy ), &
& 2._wp*ABS( zdyt (ji,jj-1) ), &
& 2._wp*ABS( zdyt (ji,jj ) ) )
END_2D
!!gm + !!st ticket ? comparaison pommes et carrottes ABS(zzslpx) et 2._wp*ABS( zdxt (ji-1,jj) )
!
DO_2D( 1, 0, 1, 0 ) !-- MUSCL horizontal advective fluxes
z0u = SIGN( 0.5_wp, pU(ji,jj,jk) )
zalpha = 0.5_wp - z0u
zu = z0u - 0.5_wp * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj)
zzwy = pt(ji ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji ,jj)
zwx(ji,jj) = pU(ji,jj,jk) * ( zalpha * zzwx + (1._wp-zalpha) * zzwy )
!
z0v = SIGN( 0.5_wp, pV(ji,jj,jk) )
zalpha = 0.5_wp - z0v
zv = z0v - 0.5_wp * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1)
zzwy = pt(ji,jj ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj )
zwy(ji,jj) = pV(ji,jj,jk) * ( zalpha * zzwx + (1._wp-zalpha) * zzwy )
END_2D
!
DO_2D( 0, 0, 0, 0 ) !-- Tracer advective trend
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj) - zwx(ji-1,jj ) &
& + zwy(ji,jj) - zwy(ji ,jj-1) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_2D
END DO
!!gm + !!st to be done with the whole rewritting of trd
!! trd routine arguments MUST be changed adding jk and zwx, zwy in 2D
!!
!! ! ! trend diagnostics
!! IF( l_trd ) THEN
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jk, jptra_xad, zwx(:,:), pU, pt(:,:,:,jn,Kbb) )
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jk, jptra_yad, zwy(:,:), pV, pt(:,:,:,jn,Kbb) )
!! END IF
!! ! ! "Poleward" heat and salt transports
!! IF( l_ptr ) CALL dia_ptr_hst( jn, jk, 'adv', zwy(:,:) )
!! ! ! heat transport
!! IF( l_hst ) CALL dia_ar5_hst( jn, jk, 'adv', zwx(:,:), zwy(:,:) )
!
! !* Vertical advective fluxes
!
#define zdzt_kp1 zdxt
#define zslpz zslpx
#define zfW zwx
zfW (A2D(0)) = 0._wp ! anciennement zwx at jk = 1
! ! anciennement zwx at jk = 2
DO_2D( 0, 0, 0, 0 )
zdzt_kp1(ji,jj) = tmask(ji,jj,2) * ( pt(ji,jj,1,jn,Kbb) - pt(ji,jj,2,jn,Kbb) )
END_2D
zslpz (A2D(0)) = 0._wp ! anciennement zslpx at jk = 1
!
IF( ln_linssh ) THEN !-- linear ssh : non zero top values
DO_2D( 0, 0, 0, 0 ) ! at the ocean surface
zfW(ji,jj) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) ! surface flux
END_2D
IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean)
DO_2D( 0, 0, 0, 0 ) ! update pt(Krhs) under the ice-shelf
ik = mikt(ji,jj) ! the flux at ik-1 is zero ( inside ice-shelf )
IF( ik > 1 ) THEN
pt(ji,jj,ik,jn,Krhs) = pt(ji,jj,ik,jn,Krhs) - pW(ji,jj,ik) * pt(ji,jj,ik,jn,Kbb) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
ENDIF
END_2D
ENDIF
ENDIF
!
! wmask usage for computing zw and zwk is needed in isf case and linear ssh
!
!
DO jk = 1, jpkm1
IF( jk < jpkm1 ) THEN
DO_2D( 0, 0, 0, 0 )
! !-- Slopes of tracer
! ! masked vertical gradient at jk+2
zdzt_kp2 = ( pt(ji,jj,jk+1,jn,Kbb) - pt(ji,jj,jk+2,jn,Kbb) ) * tmask(ji,jj,jk+2) !!st wmask(ji,jj,jk+2)
! ! vertical slope at jk+1
zslpz_kp1 = ( zdzt_kp1(ji,jj) + zdzt_kp2 ) &
& * ( 0.25_wp + SIGN( 0.25_wp, zdzt_kp1(ji,jj) * zdzt_kp2 ) )
! ! slope limitation
zslpz_kp1 = SIGN( 1.0_wp, zslpz_kp1 ) * MIN( ABS( zslpz_kp1 ), &
& 2.*ABS( zdzt_kp2 ), &
& 2.*ABS( zdzt_kp1(ji,jj) ) )
! !-- vertical advective flux at jk+1
! ! caution: zfW_kp1 is masked for ice-shelf cavities
! ! since top fluxes already added to pt(Krhs) before the vertical loop
z0w = SIGN( 0.5_wp, pW(ji,jj,jk+1) )
zalpha = 0.5_wp + z0w
zw = z0w - 0.5_wp * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpz_kp1
zzwy = pt(ji,jj,jk ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpz(ji,jj)
zfW_kp1 = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)!!st * wmask(ji,jj,jk+1)
! !-- vertical advective trend at jk
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zfW(ji,jj) - zfW_kp1 ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
! ! updates for next level
zdzt_kp1(ji,jj) = zdzt_kp2
zslpz (ji,jj) = zslpz_kp1
zfW (ji,jj) = zfW_kp1
END_2D
ELSE
DO_2D( 0, 0, 0, 0 ) !-- vertical advective trend at jpkm1
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zfW(ji,jj) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_2D
ENDIF
END DO ! end of jk loop
!
!!gm + !!st idem see above
!! ! ! send trends for diagnostic
!! IF( l_trd ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) )
!
END DO ! end of tracer loop
!
END SUBROUTINE tra_adv_mus
SUBROUTINE tra_adv_mus_hls1( kt, kit000, cdtype, p2dt, pU, pV, pW, &
& Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_mus ***
......@@ -178,7 +392,7 @@ CONTAINS
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Tracer advective trend
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
! ! trend diagnostics
......@@ -239,7 +453,7 @@ CONTAINS
!
END DO ! end of tracer loop
!
END SUBROUTINE tra_adv_mus
END SUBROUTINE tra_adv_mus_hls1
!!======================================================================
END MODULE traadv_mus
src/OCE/TRA/trazdf.F90
View file @ ee6ce8c5
......@@ -6,6 +6,7 @@ MODULE trazdf
!! History : 1.0 ! 2005-11 (G. Madec) Original code
!! 3.0 ! 2008-01 (C. Ethe, G. Madec) merge TRC-TRA
!! 4.0 ! 2017-06 (G. Madec) remove explict time-stepping option
!! 4.5 ! 2022-06 (G. Madec) refactoring to reduce memory usage (j-k-i loops)
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -22,7 +23,7 @@ MODULE trazdf
USE ldfslp ! lateral diffusion: iso-neutral slope
USE trd_oce ! trends: ocean variables
USE trdtra ! trends: tracer trend manager
USE eosbn2, ONLY: ln_SEOS, rn_b0
USE eosbn2 , ONLY: ln_SEOS, rn_b0
!
USE in_out_manager ! I/O manager
USE prtctl ! Print control
......@@ -77,7 +78,7 @@ CONTAINS
ENDIF
!
! !* compute lateral mixing trend and add it to the general trend
CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts )
CALL tra_zdf_imp( 'TRA', Kbb, Kmm, Krhs, pts, Kaa, jpts )
!!gm WHY here ! and I don't like that !
! DRAKKAR SSS control {
......@@ -116,7 +117,7 @@ CONTAINS
END SUBROUTINE tra_zdf
SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt )
SUBROUTINE tra_zdf_imp( cdtype, Kbb, Kmm, Krhs, pt, Kaa, kjpt )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_zdf_imp ***
!!
......@@ -136,128 +137,177 @@ CONTAINS
!!
!! ** Action : - pt(:,:,:,:,Kaa) becomes the after tracer
!!---------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices
INTEGER , INTENT(in ) :: kit000 ! first time step index
CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
INTEGER , INTENT(in ) :: kjpt ! number of tracers
REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp) :: zrhs, zzwi, zzws ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwt, zwd, zws
REAL(wp), DIMENSION(A1Di(0),jpk) :: zwi, zwt, zwd, zws
!!---------------------------------------------------------------------
!
! ! ============= !
DO jn = 1, kjpt ! tracer loop !
! ! ============= !
! Matrix construction
! --------------------
! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
!
IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR. &
& ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
! ! ================= !
DO_1Dj( 0, 0 ) ! i-k slices loop !
! ! ================= !
DO jn = 1, kjpt ! tracer loop !
! ! ================= !
!
! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN
DO_3D( 1, 1, 1, 1, 2, jpk )
zwt(ji,jj,jk) = avt(ji,jj,jk)
END_3D
ELSE
DO_3D( 1, 1, 1, 1, 2, jpk )
zwt(ji,jj,jk) = avs(ji,jj,jk)
END_3D
ENDIF
zwt(:,:,1) = 0._wp
! Matrix construction
! --------------------
! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
!
IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution
IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)
END_3D
ELSE ! standard or triad iso-neutral operator
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
END_3D
IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR. &
& ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
!
! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
!
IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ! use avt for temperature
!
IF( l_ldfslp ) THEN ! use avt + isoneutral diffusion contribution
IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avt(ji,jj,jk) + akz(ji,jj,jk)
END_2D
ELSE ! standard or triad iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
END_2D
ENDIF
ELSE ! use avt only
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avt(ji,jj,jk)
END_2D
ENDIF
!
ELSE ! use avs for salinty or passive tracers
!
IF( l_ldfslp ) THEN ! use avs + isoneutral diffusion contribution
IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avs(ji,jj,jk) + akz(ji,jj,jk)
END_2D
ELSE ! standard or triad iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avs(ji,jj,jk) + ah_wslp2(ji,jj,jk)
END_2D
ENDIF
ELSE !
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avs(ji,jj,jk)
END_2D
ENDIF
ENDIF
zwt(:,1) = 0._wp
!
! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zzwi = - rDt * zwt(ji,jk ) / e3w(ji,jj,jk ,Kmm)
zzws = - rDt * zwt(ji,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws &
& + rDt * ( MAX( wi(ji,jj,jk ) , 0._wp ) &
& - MIN( wi(ji,jj,jk+1) , 0._wp ) )
zwi(ji,jk) = zzwi + rDt * MIN( wi(ji,jj,jk ) , 0._wp )
zws(ji,jk) = zzws - rDt * MAX( wi(ji,jj,jk+1) , 0._wp )
END_2D
ELSE
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zwi(ji,jk) = - rDt * zwt(ji,jk ) / e3w(ji,jj,jk,Kmm)
zws(ji,jk) = - rDt * zwt(ji,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jk) - zws(ji,jk)
END_2D
ENDIF
!
!!gm BUG?? : if edmfm is equivalent to a w ==>>> just add +/- rDt * edmfm(ji,jj,jk+1/jk )
!! but edmfm is at t-point !!!! crazy??? why not keep it at w-point????
!
IF( ln_zdfmfc ) THEN ! add upward Mass Flux in the matrix
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zws(ji,jk) = zws(ji,jk) + e3t(ji,jj,jk,Kaa) * rDt * edmfm(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jk) = zwd(ji,jk) - e3t(ji,jj,jk,Kaa) * rDt * edmfm(ji,jj,jk ) / e3w(ji,jj,jk+1,Kmm)
END_2D
ENDIF
! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
! edmfa(ji,jj,jk) = 0._wp
! edmfb(ji,jj,jk) = -edmfm(ji,jj,jk ) / e3w(ji,jj,jk+1,Kmm)
! edmfc(ji,jj,jk) = edmfm(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
! END_3D
!!gm BUG : level jpk never used in the inversion
! DO_2D( 0, 0, 0, 0 )
! edmfa(ji,jj,jpk) = -edmfm(ji,jj,jpk-1) / e3w(ji,jj,jpk,Kmm)
! edmfb(ji,jj,jpk) = edmfm(ji,jj,jpk ) / e3w(ji,jj,jpk,Kmm)
! edmfc(ji,jj,jpk) = 0._wp
! END_2D
!!
!!gm BUG ??? below e3t_Kmm should be used ?
!! or even no multiplication by e3t unless there is a bug in wi calculation
!!
! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
!!gm edmfa = 0._wp except at jpk which is not used ==>> zdiagi update is useless !
! zdiagi(ji,jj,jk) = zdiagi(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfa(ji,jj,jk)
! zdiags(ji,jj,jk) = zdiags(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfc(ji,jj,jk)
! zdiagd(ji,jj,jk) = zdiagd(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfb(ji,jj,jk)
! END_3D
!!gm CALL diag_mfc( zwi, zwd, zws, rDt, Kaa )
!!gm SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa )
!
!! Matrix inversion from the first level
!!----------------------------------------------------------------------
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
!
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ( ... )( ... ) ( ... )
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
!
! m is decomposed in the product of an upper and lower triangular matrix.
! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
! and "superior" (above diagonal) components of the tridiagonal system.
! The solution will be in the 4d array pta.
! The 3d array zwt is used as a work space array.
! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then
! used as a work space array: its value is modified.
!
DO_1Di( 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) ! done one for all passive tracers (so included in the IF instruction)
zwt(ji,1) = zwd(ji,1)
END_1D
DO_2Dik( 0, 0, 2, jpkm1, 1 )
zwt(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwt(ji,jk-1)
END_2D
!
ENDIF
!
! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm)
zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws &
& + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp )
zws(ji,jj,jk) = zzws - p2dt * MAX( wi(ji,jj,jk+1) , 0._wp )
END_3D
ELSE
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm)
zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk)
END_3D
IF( ln_zdfmfc ) THEN ! add Mass Flux to the RHS
DO_2Dik( 0, 0, 1, jpkm1, 1 )
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + edmftra(ji,jj,jk,jn)
END_2D
!!gm CALL rhs_mfc( pt(:,:,:,jn,Krhs), jn )
ENDIF
!
! Modification of diagonal to add MF scheme
IF ( ln_zdfmfc ) THEN
CALL diag_mfc( zwi, zwd, zws, p2dt, Kaa )
END IF
!
!! Matrix inversion from the first level
!!----------------------------------------------------------------------
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
!
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ( ... )( ... ) ( ... )
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
!
! m is decomposed in the product of an upper and lower triangular matrix.
! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
! and "superior" (above diagonal) components of the tridiagonal system.
! The solution will be in the 4d array pta.
! The 3d array zwt is used as a work space array.
! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then
! used as a work space array: its value is modified.
!
DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) ! done one for all passive tracers (so included in the IF instruction)
zwt(ji,jj,1) = zwd(ji,jj,1)
DO_1Di( 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb ) &
& + rDt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
END_1D
DO_2Dik( 0, 0, 2, jpkm1, 1 )
zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb ) &
& + rDt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side
pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jk) / zwt(ji,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
END_2D
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
END_3D
!
ENDIF
!
! Modification of rhs to add MF scheme
IF ( ln_zdfmfc ) THEN
CALL rhs_mfc( pt(:,:,:,jn,Krhs), jn )
END IF
!
DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) &
& + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
END_2D
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) &
& + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side
pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
END_3D
!
DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer)
pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
END_2D
DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) &
& / zwt(ji,jj,jk) * tmask(ji,jj,jk)
END_3D
DO_1Di( 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer)
pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jpkm1) * tmask(ji,jj,jpkm1)
END_1D
DO_2Dik( 0, 0, jpk-2, 1, -1 )
pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jk) * pt(ji,jj,jk+1,jn,Kaa) ) &
& / zwt(ji,jk) * tmask(ji,jj,jk)
END_2D
! ! ================= !
END DO ! tracer loop !
! ! ================= !
END DO ! end tracer loop !
END_1D ! i-k slices loop !
! ! ================= !
END SUBROUTINE tra_zdf_imp
......
src/OCE/TRD/trddyn.F90
View file @ ee6ce8c5
......@@ -56,10 +56,13 @@ CONTAINS
INTEGER , INTENT(in ) :: ktrd ! trend index
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! time level index
INTEGER :: ji, jj, jk ! lopp indices
!!----------------------------------------------------------------------
!
putrd(:,:,:) = putrd(:,:,:) * umask(:,:,:) ! mask the trends
pvtrd(:,:,:) = pvtrd(:,:,:) * vmask(:,:,:)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
putrd(ji,jj,jk) = putrd(ji,jj,jk) * umask(ji,jj,jk) ! mask the trends
pvtrd(ji,jj,jk) = pvtrd(ji,jj,jk) * vmask(ji,jj,jk)
END_3D
!
!!gm NB : here a lbc_lnk should probably be added
......@@ -120,10 +123,10 @@ CONTAINS
CALL iom_put( "vtrd_rvo", pvtrd )
CASE( jpdyn_keg ) ; CALL iom_put( "utrd_keg", putrd ) ! Kinetic Energy gradient (or had)
CALL iom_put( "vtrd_keg", pvtrd )
ALLOCATE( z3dx(jpi,jpj,jpk) , z3dy(jpi,jpj,jpk) )
z3dx(:,:,:) = 0._wp ! U.dxU & V.dyV (approximation)
z3dy(:,:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! no mask as un,vn are masked
ALLOCATE( z3dx(A2D(0),jpk) , z3dy(A2D(0),jpk) ) ! U.dxU & V.dyV (approximation)
z3dx(A2D(0),jpk) = 0._wp
z3dy(A2D(0),jpk) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! no mask as un,vn are masked
z3dx(ji,jj,jk) = uu(ji,jj,jk,Kmm) * ( uu(ji+1,jj,jk,Kmm) - uu(ji-1,jj,jk,Kmm) ) / ( 2._wp * e1u(ji,jj) )
z3dy(ji,jj,jk) = vv(ji,jj,jk,Kmm) * ( vv(ji,jj+1,jk,Kmm) - vv(ji,jj-1,jk,Kmm) ) / ( 2._wp * e2v(ji,jj) )
END_3D
......@@ -139,9 +142,11 @@ CONTAINS
CALL iom_put( "vtrd_zdf", pvtrd )
!
! ! wind stress trends
ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) )
z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( e3u(:,:,1,Kmm) * rho0 )
z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( e3v(:,:,1,Kmm) * rho0 )
ALLOCATE( z2dx(A2D(0)) , z2dy(A2D(0)) )
DO_2D( 0, 0, 0, 0 )
z2dx(ji,jj) = ( utau_b(ji,jj) + utauU(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
z2dy(ji,jj) = ( vtau_b(ji,jj) + vtauV(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
END_2D
CALL iom_put( "utrd_tau", z2dx )
CALL iom_put( "vtrd_tau", z2dy )
DEALLOCATE( z2dx , z2dy )
......
src/OCE/TRD/trdglo.F90
View file @ ee6ce8c5
......@@ -42,6 +42,7 @@ MODULE trdglo
REAL(wp) :: tvolv ! volume of the whole ocean computed at v-points
REAL(wp) :: rpktrd ! potential to kinetic energy conversion
REAL(wp) :: peke ! conversion potential energy - kinetic energy trend
REAL(wp) :: xcof
! !!! domain averaged trends
REAL(wp), DIMENSION(jptot_tra) :: tmo, smo ! temperature and salinity trends
......@@ -77,44 +78,47 @@ CONTAINS
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ikbu, ikbv ! local integers
REAL(wp):: zvm, zvt, zvs, z1_2rho0 ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: ztswu, ztswv, z2dx, z2dy ! 2D workspace
REAL(wp), DIMENSION(A2D(0)) :: ztswu, ztswv, z2dx, z2dy ! 2D workspace
!!----------------------------------------------------------------------
!
IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
!
SELECT CASE( ctype )
!
CASE( 'TRA' ) !== Tracers (T & S) ==!
DO_3D( 1, 1, 1, 1, 1, jpkm1 ) ! global sum of mask volume trend and trend*T (including interior mask)
zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
zvt = ptrdx(ji,jj,jk) * zvm
zvs = ptrdy(ji,jj,jk) * zvm
tmo(ktrd) = tmo(ktrd) + zvt
smo(ktrd) = smo(ktrd) + zvs
t2 (ktrd) = t2(ktrd) + zvt * ts(ji,jj,jk,jp_tem,Kmm)
s2 (ktrd) = s2(ktrd) + zvs * ts(ji,jj,jk,jp_sal,Kmm)
END_3D
SELECT CASE( ctype )
!
CASE( 'TRA' ) !== Tracers (T & S) ==!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! global sum of mask volume trend and trend*T (including interior mask)
zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
zvt = ptrdx(ji,jj,jk) * zvm
zvs = ptrdy(ji,jj,jk) * zvm
tmo(ktrd) = tmo(ktrd) + zvt
smo(ktrd) = smo(ktrd) + zvs
t2 (ktrd) = t2(ktrd) + zvt * ts(ji,jj,jk,jp_tem,Kmm)
s2 (ktrd) = s2(ktrd) + zvs * ts(ji,jj,jk,jp_sal,Kmm)
END_3D
! ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
IF( ln_linssh .AND. ktrd == jptra_zad ) THEN
tmo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
smo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
t2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
s2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
ENDIF
IF( ln_linssh .AND. ktrd == jptra_zad ) THEN
DO_2D( 0, 0, 0, 0 ) ! global sum of mask volume trend and trend*T (including interior mask)
zvm = ww(ji,jj,1) * e1e2t(ji,jj) * tmask_i(ji,jj)
tmo(jptra_sad) = tmo(jptra_sad) + ts(ji,jj,1,jp_tem,Kmm) * zvm
smo(jptra_sad) = smo(jptra_sad) + ts(ji,jj,1,jp_sal,Kmm) * zvm
t2 (jptra_sad) = t2 (jptra_sad) + ts(ji,jj,1,jp_tem,Kmm) * ts(ji,jj,1,jp_tem,Kmm) * zvm
s2 (jptra_sad) = s2 (jptra_sad) + ts(ji,jj,1,jp_sal,Kmm) * ts(ji,jj,1,jp_sal,Kmm) * zvm
END_2D
ENDIF
!
IF( ktrd == jptra_atf ) THEN ! last trend (asselin time filter)
!
CALL glo_tra_wri( kt ) ! print the results in ocean.output
!
tmo(:) = 0._wp ! prepare the next time step (domain averaged array reset to zero)
smo(:) = 0._wp
t2 (:) = 0._wp
s2 (:) = 0._wp
!
ENDIF
IF( ktrd == jptra_atf ) THEN ! last trend (asselin time filter)
!
CALL glo_tra_wri( kt ) ! print the results in ocean.output
!
tmo(:) = 0._wp ! prepare the next time step (domain averaged array reset to zero)
smo(:) = 0._wp
t2 (:) = 0._wp
s2 (:) = 0._wp
!
ENDIF
!
CASE( 'DYN' ) !== Momentum and KE ==!
DO_3D( 1, 0, 1, 0, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk) &
& * e1e2u (ji,jj) * e3u(ji,jj,jk,Kmm)
zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk) &
......@@ -126,11 +130,11 @@ CONTAINS
!
IF( ktrd == jpdyn_zdf ) THEN ! zdf trend: compute separately the surface forcing trend
z1_2rho0 = 0.5_wp / rho0
DO_2D( 1, 0, 1, 0 )
zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk) &
& * z1_2rho0 * e1e2u(ji,jj)
zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk) &
& * z1_2rho0 * e1e2v(ji,jj)
DO_2D( 0, 0, 0, 0 )
zvt = ( utau_b(ji,jj) + utauU(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk) &
& * z1_2rho0 * e1e2u(ji,jj)
zvs = ( vtau_b(ji,jj) + vtauV(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk) &
& * z1_2rho0 * e1e2v(ji,jj)
umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
hke(jpdyn_tau) = hke(jpdyn_tau) + uu(ji,jj,1,Kmm) * zvt + vv(ji,jj,1,Kmm) * zvs
......@@ -184,8 +188,7 @@ CONTAINS
INTEGER, INTENT(in) :: Kmm ! time level index
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zcof ! local scalar
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zkx, zky, zkz, zkepe
REAL(wp), DIMENSION(A2D(0),jpk) :: zkx, zky, zkz, zkepe
!!----------------------------------------------------------------------
! I. Momentum trends
......@@ -196,24 +199,24 @@ CONTAINS
! I.1 Conversion potential energy - kinetic energy
! --------------------------------------------------
! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
zkx (:,:,:) = 0._wp
zky (:,:,:) = 0._wp
zkz (:,:,:) = 0._wp
zkepe(:,:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! loop from top to bottom
zkx (ji,jj,jk) = 0._wp
zky (ji,jj,jk) = 0._wp
zkz (ji,jj,jk) = 0._wp
zkepe(ji,jj,jk) = 0._wp
END_3D
CALL eos( ts(:,:,:,:,Kmm), rhd, rhop ) ! now potential density
zcof = 0.5_wp / rho0 ! Density flux at w-point
zkz(:,:,1) = 0._wp
DO jk = 2, jpk
zkz(:,:,jk) = zcof * e1e2t(:,:) * ww(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
END DO
zkz(A2D(0),1) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpk ) ! loop from top to bottom
zkz(ji,jj,jk) = xcof * e1e2t(ji,jj) * ww(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) ) * tmask_i(ji,jj)
END_3D
zcof = 0.5_wp / rho0 ! Density flux at u and v-points
DO_3D( 1, 0, 1, 0, 1, jpkm1 )
zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) &
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zkx(ji,jj,jk) = xcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) &
& * uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) &
zky(ji,jj,jk) = xcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) &
& * vv(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
END_3D
......@@ -227,10 +230,9 @@ CONTAINS
! I.2 Basin averaged kinetic energy trend
! ----------------------------------------
peke = 0._wp
DO jk = 1, jpkm1
peke = peke + SUM( zkepe(:,:,jk) * gdept(:,:,jk,Kmm) * e1e2t(:,:) &
& * e3t(:,:,jk,Kmm) )
END DO
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! loop from top to bottom
peke = peke + zkepe(ji,jj,jk) * gdept(ji,jj,jk,Kmm) * e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)
END_3D
peke = grav * peke
! I.3 Sums over the global domain
......@@ -509,11 +511,14 @@ CONTAINS
WRITE(numout,*) '~~~~~~~~~~~~~'
ENDIF
xcof = 0.5_wp / rho0
! Total volume at t-points:
tvolt = 0._wp
DO jk = 1, jpkm1
tvolt = tvolt + SUM( e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) * tmask_i(:,:) )
END DO
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Density flux divergence at t-point
tvolt = tvolt + e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
END_3D
CALL mpp_sum( 'trdglo', tvolt ) ! sum over the global domain
IF(lwp) WRITE(numout,*) ' total ocean volume at T-point tvolt = ',tvolt
......
src/OCE/TRD/trdken.F90
View file @ ee6ce8c5
......@@ -30,6 +30,10 @@ MODULE trdken
IMPLICIT NONE
PRIVATE
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
PUBLIC trd_ken ! called by trddyn module
PUBLIC trd_ken_init ! called by trdini module
......@@ -38,9 +42,6 @@ MODULE trdken
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: bu, bv ! volume of u- and v-boxes
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: r1_bt ! inverse of t-box volume
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: trdken.F90 15104 2021-07-07 14:36:00Z clem $
......@@ -52,7 +53,7 @@ CONTAINS
!!---------------------------------------------------------------------
!! *** FUNCTION trd_ken_alloc ***
!!---------------------------------------------------------------------
ALLOCATE( bu(jpi,jpj,jpk) , bv(jpi,jpj,jpk) , r1_bt(jpi,jpj,jpk) , STAT= trd_ken_alloc )
ALLOCATE( bu(A2D(0),jpk) , bv(A2D(0),jpk) , r1_bt(A2D(0),jpk) , STAT= trd_ken_alloc )
!
CALL mpp_sum ( 'trdken', trd_ken_alloc )
IF( trd_ken_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trd_ken_alloc: failed to allocate arrays' )
......@@ -77,31 +78,32 @@ CONTAINS
!
!
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: putrd, pvtrd ! U and V masked trends
INTEGER , INTENT(in ) :: ktrd ! trend index
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: putrd, pvtrd ! U and V masked trends
INTEGER , INTENT(in ) :: ktrd ! trend index
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! time level index
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ikbu , ikbv ! local integers
INTEGER :: ikbum1, ikbvm1 ! - -
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z2dx, z2dy, zke2d ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zke ! 3D workspace
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zke2d ! 2D workspace
REAL(wp) :: z2dx, z2dy, z2dxm1, z2dym1 ! 2D workspace
REAL(wp), DIMENSION(A2D(0),jpk) :: zke ! 3D workspace
!!----------------------------------------------------------------------
!
CALL lbc_lnk( 'trdken', putrd, 'U', -1.0_wp , pvtrd, 'V', -1.0_wp ) ! lateral boundary conditions
!
nkstp = kt
DO jk = 1, jpkm1
bu (:,:,jk) = e1e2u(:,:) * e3u(:,:,jk,Kmm)
bv (:,:,jk) = e1e2v(:,:) * e3v(:,:,jk,Kmm)
r1_bt(:,:,jk) = r1_e1e2t(:,:) / e3t(:,:,jk,Kmm) * tmask(:,:,jk)
DO_2D( 0, 0, 0, 0 )
bu (ji,jj,jk) = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
bv (ji,jj,jk) = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
r1_bt(ji,jj,jk) = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
END_2D
END DO
!
zke(:,:,jpk) = 0._wp
zke(1:nn_hls,:, : ) = 0._wp
zke(:,1:nn_hls, : ) = 0._wp
DO_3D( 0, nn_hls, 0, nn_hls, 1, jpkm1 )
zke(A2D(0),:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zke(ji,jj,jk) = 0.5_wp * rho0 *( uu(ji ,jj,jk,Kmm) * putrd(ji ,jj,jk) * bu(ji ,jj,jk) &
& + uu(ji-1,jj,jk,Kmm) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk) &
& + vv(ji,jj ,jk,Kmm) * pvtrd(ji,jj ,jk) * bv(ji,jj ,jk) &
......@@ -118,35 +120,31 @@ CONTAINS
CASE( jpdyn_ldf ) ; CALL iom_put( "ketrd_ldf" , zke ) ! lateral diffusion
CASE( jpdyn_zdf ) ; CALL iom_put( "ketrd_zdf" , zke ) ! vertical diffusion
! ! ! wind stress trends
ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) , zke2d(jpi,jpj) )
z2dx(:,:) = uu(:,:,1,Kmm) * ( utau_b(:,:) + utau(:,:) ) * e1e2u(:,:) * umask(:,:,1)
z2dy(:,:) = vv(:,:,1,Kmm) * ( vtau_b(:,:) + vtau(:,:) ) * e1e2v(:,:) * vmask(:,:,1)
zke2d(1:nn_hls,:) = 0._wp ; zke2d(:,1:nn_hls) = 0._wp
DO_2D( 0, nn_hls, 0, nn_hls )
zke2d(ji,jj) = r1_rho0 * 0.5_wp * ( z2dx(ji,jj) + z2dx(ji-1,jj) &
& + z2dy(ji,jj) + z2dy(ji,jj-1) ) * r1_bt(ji,jj,1)
END_2D
ALLOCATE( zke2d(A2D(0)) ) ; zke2d(A2D(0)) = 0._wp
DO_2D( 0, 0, 0, 0 )
z2dx = uu(ji,jj,1,Kmm) * ( utau_b(ji,jj) + utauU(ji,jj) ) * e1e2u(ji,jj) * umask(ji,jj,1)
z2dy = vv(ji,jj,1,Kmm) * ( vtau_b(ji,jj) + vtauV(ji,jj) ) * e1e2v(ji,jj) * vmask(ji,jj,1)
z2dxm1 = uu(ji-1,jj,1,Kmm) * ( utau_b(ji-1,jj) + utauU(ji-1,jj) ) * e1e2u(ji-1,jj) * umask(ji-1,jj,1)
z2dym1 = vv(ji,jj-1,1,Kmm) * ( vtau_b(ji,jj-1) + vtauV(ji,jj-1) ) * e1e2v(ji,jj-1) * vmask(ji,jj-1,1)
zke2d(ji,jj) = r1_rho0 * 0.5_wp * ( z2dx + z2dxm1 + z2dy + z2dym1 ) * r1_bt(ji,jj,1)
END_2D
CALL iom_put( "ketrd_tau" , zke2d ) !
DEALLOCATE( z2dx , z2dy , zke2d )
DEALLOCATE( zke2d )
CASE( jpdyn_bfr ) ; CALL iom_put( "ketrd_bfr" , zke ) ! bottom friction (explicit case)
!!gm TO BE DONE properly
!!gm only valid if ln_drgimp=F otherwise the bottom stress as to be recomputed at the end of the computation....
! IF(.NOT. ln_drgimp) THEN
! DO jj = 1, jpj !
! DO ji = 1, jpi
! ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels
! ikbv = mbkv(ji,jj)
! z2dx(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm)
! z2dy(ji,jj) = vv(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm)
! END DO
! END DO
! zke2d(A2D(0)) = 0._wp
! DO_2D( 0, 0, 0, 0 )
! ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels
! ikbv = mbkv(ji,jj)
! z2dx = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm)
! z2dy = vv(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm)
! z2dxm1 = uu(ji-1,jj,ikbu,Kmm) * bfrua(ji-1,jj) * uu(ji-1,jj,ikbu,Kmm)
! z2dym1 = vv(ji,jj-1,ikbu,Kmm) * bfrva(ji,jj-1) * vv(ji,jj-1,ikbv,Kmm)
! zke2d(ji,jj) = 0.5_wp * ( z2dx + z2dxm1 + z2dy + z2dym1 ) * r1_bt(ji,jj,1)
! END_2D
! zke2d(1,:) = 0._wp ; zke2d(:,1) = 0._wp
! DO jj = 2, jpj
! DO ji = 2, jpi
! zke2d(ji,jj) = 0.5_wp * ( z2dx(ji,jj) + z2dx(ji-1,jj) &
! & + z2dy(ji,jj) + z2dy(ji,jj-1) ) * r1_bt(ji,jj, BEURK!!!
! END DO
! END DO
! CALL iom_put( "ketrd_bfr" , zke2d ) ! bottom friction (explicit case)
! ENDIF
!!gm end
......@@ -175,11 +173,11 @@ CONTAINS
CASE( jpdyn_ken ) ; ! kinetic energy
! called in dynnxt.F90 before asselin time filter
! with putrd=uu(Krhs) and pvtrd=vv(Krhs)
zke(:,:,:) = 0.5_wp * zke(:,:,:)
zke(A2D(0),:) = 0.5_wp * zke(A2D(0),:)
CALL iom_put( "KE", zke )
!
CALL ken_p2k( kt , zke, Kmm )
CALL iom_put( "ketrd_convP2K", zke ) ! conversion -rau*g*w
CALL iom_put( "ketrd_convP2K", zke ) ! conversion -rau*g*w
!
END SELECT
!
......@@ -203,22 +201,23 @@ CONTAINS
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: iku, ikv ! local integers
REAL(wp) :: zcoef ! local scalars
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zconv ! 3D workspace
REAL(wp), DIMENSION(A2D(0),jpk) :: zconv ! 3D workspace
!!----------------------------------------------------------------------
!
! Local constant initialization
zcoef = - rho0 * grav * 0.5_wp
! Surface value (also valid in partial step case)
zconv(:,:,1) = zcoef * ( 2._wp * rhd(:,:,1) ) * ww(:,:,1) * e3w(:,:,1,Kmm)
DO_2D( 0, 0, 0, 0 )
zconv(ji,jj,1) = zcoef * ( 2._wp * rhd(ji,jj,1) ) * ww(ji,jj,1) * e3w(ji,jj,1,Kmm)
END_2D
! interior value (2=<jk=<jpkm1)
DO jk = 2, jpk
zconv(:,:,jk) = zcoef * ( rhd(:,:,jk) + rhd(:,:,jk-1) ) * ww(:,:,jk) * e3w(:,:,jk,Kmm)
END DO
DO_3D( 0, 0, 0, 0 , 2, jpk )
zconv(ji,jj,jk) = zcoef * ( rhd(ji,jj,jk) + rhd(ji,jj,jk-1) ) * ww(ji,jj,jk) * e3w(ji,jj,jk,Kmm)
END_3D
! conv value on T-point
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zcoef = 0.5_wp / e3t(ji,jj,jk,Kmm)
pconv(ji,jj,jk) = zcoef * ( zconv(ji,jj,jk) + zconv(ji,jj,jk+1) ) * tmask(ji,jj,jk)
END_3D
......
src/OCE/TRD/trdmxl.F90
View file @ ee6ce8c5
......@@ -45,20 +45,7 @@ MODULE trdmxl
PUBLIC trd_mxl_init ! routine called by opa.F90
PUBLIC trd_mxl_zint ! routine called by tracers routines
INTEGER :: nkstp ! current time step
!!gm to be moved from trdmxl_oce
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: hml ! ML depth (sum of e3t over nmln-1 levels) [m]
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tml , sml ! now ML averaged T & S
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tmlb_nf, smlb_nf ! not filtered before ML averaged T & S
!
!
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: hmlb, hmln ! before, now, and after Mixed Layer depths [m]
!
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tb_mlb, tb_mln ! before (not filtered) tracer averaged over before and now ML
!
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tn_mln ! now tracer averaged over now ML
!!gm end
INTEGER :: nkstp ! current time step
CHARACTER (LEN=40) :: clhstnam ! name of the trends NetCDF file
INTEGER :: nh_t, nmoymltrd
......@@ -111,44 +98,41 @@ CONTAINS
IF ( kt /= nkstp ) THEN != 1st call at kt time step =!
! !==============================!
nkstp = kt
! !== reset trend arrays to zero ==!
tmltrd(:,:,:) = 0._wp ; smltrd(:,:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpktrd )
tmltrd(ji,jj,jk) = 0._wp
smltrd(ji,jj,jk) = 0._wp
END_3D
!
wkx(:,:,:) = 0._wp !== now ML weights for vertical averaging ==!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpktrd ) ! initialize wkx with vertical scale factor in mixed-layer
! !== now ML weights for vertical averaging ==!
DO_3D( 0, 0, 0, 0, 1, jpktrd ) ! initialize wkx with vertical scale factor in mixed-layer
IF( jk - kmxln(ji,jj) < 0 ) THEN
wkx(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
ELSE
wkx(ji,jj,jk) = 0._wp
ENDIF
END_3D
!
hmxl(:,:) = 0._wp ! NOW mixed-layer depth
DO jk = 1, jpktrd
hmxl(:,:) = hmxl(:,:) + wkx(:,:,jk)
END DO
DO jk = 1, jpktrd ! integration weights
wkx(:,:,jk) = wkx(:,:,jk) / MAX( 1.e-20_wp, hmxl(:,:) ) * tmask(:,:,1)
END DO
DO_3D( 0, 0, 0, 0, 1, jpktrd )
hmxl(ji,jj) = hmxl(ji,jj) + wkx(ji,jj,jk)
wkx (ji,jj,jk) = wkx (ji,jj,jk) / MAX( 1.e-20_wp, hmxl(ji,jj) ) * tmask(ji,jj,1)
END_3D
!
! !== Vertically averaged T and S ==!
tml(:,:) = 0._wp ; sml(:,:) = 0._wp
DO jk = 1, jpktrd
tml(:,:) = tml(:,:) + wkx(:,:,jk) * ts(:,:,jk,jp_tem,Kmm)
sml(:,:) = sml(:,:) + wkx(:,:,jk) * ts(:,:,jk,jp_sal,Kmm)
END DO
DO_3D( 0, 0, 0, 0, 1, jpktrd )
tml(ji,jj) = tml(ji,jj) + wkx(ji,jj,jk) * ts(ji,jj,jk,jp_tem,Kmm)
sml(ji,jj) = sml(ji,jj) + wkx(ji,jj,jk) * ts(ji,jj,jk,jp_sal,Kmm)
END_3D
!
ENDIF
! mean now trends over the now ML
tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + ptrdx(:,:,jk) * wkx(:,:,jk) ! temperature
smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + ptrdy(:,:,jk) * wkx(:,:,jk) ! salinity
DO_2D( 0, 0, 0, 0 )
tmltrd(ji,jj,ktrd) = tmltrd(ji,jj,ktrd) + ptrdx(ji,jj,jk) * wkx(ji,jj,jk) ! temperature
smltrd(ji,jj,ktrd) = smltrd(ji,jj,ktrd) + ptrdy(ji,jj,jk) * wkx(ji,jj,jk) ! salinity
END_2D
!!gm to be put juste before the output !
......@@ -249,11 +233,11 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: ktrd ! ocean trend index
CHARACTER(len=2) , INTENT( in ) :: ctype ! 2D surface/bottom or 3D interior physics
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pttrdmxl ! temperature trend
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pstrdmxl ! salinity trend
REAL(wp), DIMENSION(A2D(0),jpk), INTENT( in ) :: pttrdmxl ! temperature trend
REAL(wp), DIMENSION(A2D(0),jpk), INTENT( in ) :: pstrdmxl ! salinity trend
!
INTEGER :: ji, jj, jk, isum
REAL(wp), DIMENSION(jpi,jpj) :: zvlmsk
REAL(wp), DIMENSION(A2D(0)) :: zvlmsk
!!----------------------------------------------------------------------
! I. Definition of control surface and associated fields
......@@ -268,11 +252,13 @@ CONTAINS
! ... Set nmxl(ji,jj) = index of first T point below control surf. or outside mixed-layer
IF( nn_ctls == 0 ) THEN ! * control surface = mixed-layer with density criterion
nmxl(:,:) = nmln(:,:) ! array nmln computed in zdfmxl.F90
ELSEIF( nn_ctls == 1 ) THEN ! * control surface = read index from file
IF( nn_ctls == 0 ) THEN ! * control surface = mixed-layer with density criterion
DO_2D( 0, 0, 0, 0 )
nmxl(ji,jj) = nmln(ji,jj) ! array nmln computed in zdfmxl.F90
END_2D
ELSEIF( nn_ctls == 1 ) THEN ! * control surface = read index from file
nmxl(:,:) = nbol(:,:)
ELSEIF( nn_ctls >= 2 ) THEN ! * control surface = model level
ELSEIF( nn_ctls >= 2 ) THEN ! * control surface = model level
nn_ctls = MIN( nn_ctls, jpktrd - 1 )
nmxl(:,:) = nn_ctls + 1
ENDIF
......@@ -335,9 +321,9 @@ CONTAINS
REAL(wp) :: zavt, zfn, zfn2
! ! z(ts)mltot : dT/dt over the anlysis window (including Asselin)
! ! z(ts)mlres : residual = dh/dt entrainment term
REAL(wp), DIMENSION(jpi,jpj ) :: ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf
REAL(wp), DIMENSION(jpi,jpj ) :: ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2
REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmltrd2, zsmltrd2 ! only needed for mean diagnostics
REAL(wp), DIMENSION(A2D(0) ) :: ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf
REAL(wp), DIMENSION(A2D(0) ) :: ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2
REAL(wp), DIMENSION(A2D(0),jpk) :: ztmltrd2, zsmltrd2 ! only needed for mean diagnostics
!!----------------------------------------------------------------------
! ======================================================================
......@@ -446,12 +432,7 @@ CONTAINS
itmod = kt - nit000 + 1
MODULO_NTRD : IF( MOD( itmod, nn_trd ) == 0 ) THEN ! nitend MUST be multiple of nn_trd
!
ztmltot (:,:) = 0.e0 ; zsmltot (:,:) = 0.e0 ! reset arrays to zero
ztmlres (:,:) = 0.e0 ; zsmlres (:,:) = 0.e0
ztmltot2(:,:) = 0.e0 ; zsmltot2(:,:) = 0.e0
ztmlres2(:,:) = 0.e0 ; zsmlres2(:,:) = 0.e0
!
zfn = REAL( nmoymltrd, wp ) ; zfn2 = zfn * zfn
! III.1 Prepare fields for output ("instantaneous" diagnostics)
......@@ -469,25 +450,18 @@ CONTAINS
ztmlatf(:,:) = tmlatfm(:,:) - tmlatfn(:,:) + tmlatfb(:,:)
zsmlatf(:,:) = smlatfm(:,:) - smlatfn(:,:) + smlatfb(:,:)
!-- Lateral boundary conditions
! ... temperature ... ... salinity ...
CALL lbc_lnk( 'trdmxl', ztmltot , 'T', 1.0_wp, zsmltot , 'T', 1.0_wp, &
& ztmlres , 'T', 1.0_wp, zsmlres , 'T', 1.0_wp, &
& ztmlatf , 'T', 1.0_wp, zsmlatf , 'T', 1.0_wp )
! III.2 Prepare fields for output ("mean" diagnostics)
! ----------------------------------------------------
!-- Update the ML depth time sum (to build the Leap-Frog time mean)
hmxl_sum(:,:) = hmxlbn(:,:) + 2 * ( hmxl_sum(:,:) - hmxl(:,:) ) + hmxl(:,:)
hmxl_sum(:,:) = hmxlbn(:,:) + 2._wp * ( hmxl_sum(:,:) - hmxl(:,:) ) + hmxl(:,:)
!-- Compute temperature total trends
tml_sum (:,:) = tmlbn(:,:) + 2 * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:)
tml_sum (:,:) = tmlbn(:,:) + 2._wp * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:)
ztmltot2(:,:) = ( tml_sum(:,:) - tml_sumb(:,:) ) / p2dt ! now in degC/s
!-- Compute salinity total trends
sml_sum (:,:) = smlbn(:,:) + 2 * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:)
sml_sum (:,:) = smlbn(:,:) + 2._wp * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:)
zsmltot2(:,:) = ( sml_sum(:,:) - sml_sumb(:,:) ) / p2dt ! now in psu/s
!-- Compute temperature residuals
......@@ -495,7 +469,7 @@ CONTAINS
ztmltrd2(:,:,jl) = tmltrd_csum_ub(:,:,jl) + tmltrd_csum_ln(:,:,jl)
END DO
ztmltrdm2(:,:) = 0.e0
ztmltrdm2(:,:) = 0._wp
DO jl = 1, jpltrd
ztmltrdm2(:,:) = ztmltrdm2(:,:) + ztmltrd2(:,:,jl)
END DO
......@@ -508,7 +482,7 @@ CONTAINS
zsmltrd2(:,:,jl) = smltrd_csum_ub(:,:,jl) + smltrd_csum_ln(:,:,jl)
END DO
zsmltrdm2(:,:) = 0.
zsmltrdm2(:,:) = 0._wp
DO jl = 1, jpltrd
zsmltrdm2(:,:) = zsmltrdm2(:,:) + zsmltrd2(:,:,jl)
END DO
......@@ -519,13 +493,6 @@ CONTAINS
!-- Diagnose Asselin trend over the analysis window
ztmlatf2(:,:) = ztmltrd2(:,:,jpmxl_atf) - tmltrd_sum(:,:,jpmxl_atf) + tmltrd_atf_sumb(:,:)
zsmlatf2(:,:) = zsmltrd2(:,:,jpmxl_atf) - smltrd_sum(:,:,jpmxl_atf) + smltrd_atf_sumb(:,:)
!-- Lateral boundary conditions
! ... temperature ... ... salinity ...
CALL lbc_lnk( 'trdmxl', ztmltot2, 'T', 1.0_wp, zsmltot2, 'T', 1.0_wp, &
& ztmlres2, 'T', 1.0_wp, zsmlres2, 'T', 1.0_wp )
!
CALL lbc_lnk( 'trdmxl', ztmltrd2(:,:,:), 'T', 1.0_wp, zsmltrd2(:,:,:), 'T', 1.0_wp ) ! / in the NetCDF trends file
! III.3 Time evolution array swap
! -------------------------------
......@@ -622,6 +589,8 @@ CONTAINS
! ======================================================================
!-- Write the trends for T/S instantaneous diagnostics
! clem => these fields do not exist in field_def
IF( ln_trdmxl_instant ) THEN
......
src/OCE/TRD/trdmxl_oce.F90
View file @ ee6ce8c5
......@@ -80,6 +80,8 @@ MODULE trdmxl_oce
smltrd_csum_ln, & !: ( idem for salinity )
smltrd_csum_ub !:
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: trdmxl_oce.F90 10425 2018-12-19 21:54:16Z smasson $
......@@ -101,29 +103,29 @@ CONTAINS
ierr(:) = 0
ALLOCATE( nmxl (jpi,jpj) , nbol (jpi,jpj), &
& wkx (jpi,jpj,jpk), hmxl (jpi,jpj), &
& tml (jpi,jpj) , sml (jpi,jpj), &
& tmlb (jpi,jpj) , smlb (jpi,jpj), &
& tmlbb(jpi,jpj) , smlbb(jpi,jpj), STAT = ierr(1) )
ALLOCATE( tmlbn(jpi,jpj) , smlbn(jpi,jpj), &
& tmltrdm(jpi,jpj), smltrdm(jpi,jpj), &
& tml_sum(jpi,jpj), tml_sumb(jpi,jpj),&
& tmltrd_atf_sumb(jpi,jpj) , STAT=ierr(2) )
ALLOCATE( sml_sum(jpi,jpj), sml_sumb(jpi,jpj), &
& smltrd_atf_sumb(jpi,jpj), &
& hmxl_sum(jpi,jpj), hmxlbn(jpi,jpj), &
& tmlatfb(jpi,jpj), tmlatfn(jpi,jpj), STAT = ierr(3) )
ALLOCATE( smlatfb(jpi,jpj), smlatfn(jpi,jpj), &
& tmlatfm(jpi,jpj), smlatfm(jpi,jpj), &
& tmltrd(jpi,jpj,jpltrd), smltrd(jpi,jpj,jpltrd), STAT=ierr(4))
ALLOCATE( tmltrd_sum(jpi,jpj,jpltrd),tmltrd_csum_ln(jpi,jpj,jpltrd), &
& tmltrd_csum_ub(jpi,jpj,jpltrd), smltrd_sum(jpi,jpj,jpltrd), &
& smltrd_csum_ln(jpi,jpj,jpltrd), smltrd_csum_ub(jpi,jpj,jpltrd), STAT=ierr(5) )
ALLOCATE( nmxl (A2D(0)) , nbol (A2D(0)), &
& wkx (A2D(0),jpk), hmxl (A2D(0)), &
& tml (A2D(0)) , sml (A2D(0)), &
& tmlb (A2D(0)) , smlb (A2D(0)), &
& tmlbb(A2D(0)) , smlbb(A2D(0)), STAT = ierr(1) )
ALLOCATE( tmlbn(A2D(0)) , smlbn(A2D(0)), &
& tmltrdm(A2D(0)), smltrdm(A2D(0)), &
& tml_sum(A2D(0)), tml_sumb(A2D(0)),&
& tmltrd_atf_sumb(A2D(0)) , STAT=ierr(2) )
ALLOCATE( sml_sum(A2D(0)), sml_sumb(A2D(0)), &
& smltrd_atf_sumb(A2D(0)), &
& hmxl_sum(A2D(0)), hmxlbn(A2D(0)), &
& tmlatfb(A2D(0)), tmlatfn(A2D(0)), STAT = ierr(3) )
ALLOCATE( smlatfb(A2D(0)), smlatfn(A2D(0)), &
& tmlatfm(A2D(0)), smlatfm(A2D(0)), &
& tmltrd(A2D(0),jpltrd), smltrd(A2D(0),jpltrd), STAT=ierr(4))
ALLOCATE( tmltrd_sum(A2D(0),jpltrd),tmltrd_csum_ln(A2D(0),jpltrd), &
& tmltrd_csum_ub(A2D(0),jpltrd), smltrd_sum(A2D(0),jpltrd), &
& smltrd_csum_ln(A2D(0),jpltrd), smltrd_csum_ub(A2D(0),jpltrd), STAT=ierr(5) )
!
trdmxl_oce_alloc = MAXVAL( ierr )
CALL mpp_sum ( 'trdmxl_oce', trdmxl_oce_alloc )
......
src/OCE/TRD/trdpen.F90
View file @ ee6ce8c5
......@@ -35,6 +35,7 @@ MODULE trdpen
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: rab_pe ! partial derivatives of PE anomaly with respect to T and S
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -48,7 +49,7 @@ CONTAINS
!!---------------------------------------------------------------------
!! *** FUNCTION trd_tra_alloc ***
!!---------------------------------------------------------------------
ALLOCATE( rab_pe(jpi,jpj,jpk,jpts) , STAT= trd_pen_alloc )
ALLOCATE( rab_pe(A2D(0),jpk,jpts) , STAT= trd_pen_alloc )
!
CALL mpp_sum ( 'trdpen', trd_pen_alloc )
IF( trd_pen_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trd_pen_alloc: failed to allocate arrays' )
......@@ -69,37 +70,40 @@ CONTAINS
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp) , INTENT(in) :: pdt ! time step [s]
!
INTEGER :: jk ! dummy loop indices
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpe ! 3D workspace
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D workspace
REAL(wp), DIMENSION(A2D(0),jpk) :: zpe ! 3D workspace
!!----------------------------------------------------------------------
!
zpe(:,:,:) = 0._wp
zpe(A2D(0),:) = 0._wp
!
IF( kt /= nkstp ) THEN ! full eos: set partial derivatives at the 1st call of kt time step
nkstp = kt
CALL eos_pen( ts(:,:,:,:,Kmm), rab_PE, zpe, Kmm )
CALL eos_pen( ts(A2D(0),:,:,Kmm), rab_pe, zpe, Kmm )
CALL iom_put( "alphaPE", rab_pe(:,:,:,jp_tem) )
CALL iom_put( "betaPE" , rab_pe(:,:,:,jp_sal) )
CALL iom_put( "PEanom" , zpe )
ENDIF
!
zpe(:,:,jpk) = 0._wp
DO jk = 1, jpkm1
zpe(:,:,jk) = ( - ( rab_n(:,:,jk,jp_tem) + rab_pe(:,:,jk,jp_tem) ) * ptrdx(:,:,jk) &
& + ( rab_n(:,:,jk,jp_sal) + rab_pe(:,:,jk,jp_sal) ) * ptrdy(:,:,jk) )
END DO
zpe(A2D(0),jpk) = 0._wp
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zpe(ji,jj,jk) = ( - ( rab_n(ji,jj,jk,jp_tem) + rab_pe(ji,jj,jk,jp_tem) ) * ptrdx(ji,jj,jk) &
& + ( rab_n(ji,jj,jk,jp_sal) + rab_pe(ji,jj,jk,jp_sal) ) * ptrdy(ji,jj,jk) )
END_3D
SELECT CASE ( ktrd )
CASE ( jptra_xad ) ; CALL iom_put( "petrd_xad", zpe ) ! zonal advection
CASE ( jptra_yad ) ; CALL iom_put( "petrd_yad", zpe ) ! merid. advection
CASE ( jptra_zad ) ; CALL iom_put( "petrd_zad", zpe ) ! vertical advection
IF( ln_linssh ) THEN ! cst volume : adv flux through z=0 surface
ALLOCATE( z2d(jpi,jpj) )
z2d(:,:) = ww(:,:,1) * ( &
& - ( rab_n(:,:,1,jp_tem) + rab_pe(:,:,1,jp_tem) ) * ts(:,:,1,jp_tem,Kmm) &
& + ( rab_n(:,:,1,jp_sal) + rab_pe(:,:,1,jp_sal) ) * ts(:,:,1,jp_sal,Kmm) &
& ) / e3t(:,:,1,Kmm)
ALLOCATE( z2d(A2D(0)) )
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = ww(ji,jj,1) * ( &
& - ( rab_n(ji,jj,1,jp_tem) + rab_pe(ji,jj,1,jp_tem) ) * ts(ji,jj,1,jp_tem,Kmm) &
& + ( rab_n(ji,jj,1,jp_sal) + rab_pe(ji,jj,1,jp_sal) ) * ts(ji,jj,1,jp_sal,Kmm) &
& ) / e3t(ji,jj,1,Kmm)
END_2D
CALL iom_put( "petrd_sad" , z2d )
DEALLOCATE( z2d )
ENDIF
......
src/OCE/TRD/trdtra.F90
View file @ ee6ce8c5
......@@ -53,7 +53,7 @@ CONTAINS
!!---------------------------------------------------------------------
!! *** FUNCTION trd_tra_alloc ***
!!---------------------------------------------------------------------
ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , avt_evd(jpi,jpj,jpk), STAT= trd_tra_alloc )
ALLOCATE( trdtx(A2D(0),jpk), trdty(A2D(0),jpk), trdt(A2D(0),jpk), avt_evd(A2D(0),jpk), STAT= trd_tra_alloc )
!
CALL mpp_sum ( 'trdtra', trd_tra_alloc )
IF( trd_tra_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trd_tra_alloc: failed to allocate arrays' )
......@@ -73,24 +73,25 @@ CONTAINS
!! - 'TRA' case : regroup T & S trends
!! - send the trends to trd_tra_mng (trdtrc) for further processing
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! time step
CHARACTER(len=3) , INTENT(in) :: ctype ! tracers trends type 'TRA'/'TRC'
INTEGER , INTENT(in) :: ktra ! tracer index
INTEGER , INTENT(in) :: ktrd ! tracer trend index
INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptrd ! tracer trend or flux
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: pu ! now velocity
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: ptra ! now tracer variable
INTEGER , INTENT(in) :: kt ! time step
CHARACTER(len=3) , INTENT(in) :: ctype ! tracers trends type 'TRA'/'TRC'
INTEGER , INTENT(in) :: ktra ! tracer index
INTEGER , INTENT(in) :: ktrd ! tracer trend index
INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices
REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptrd ! tracer trend or flux
REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: pu ! now velocity
REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: ptra ! now tracer variable
!
INTEGER :: jk ! loop indices
INTEGER :: i01 ! 0 or 1
REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztrds ! 3D workspace
INTEGER :: ji,jj,jk ! loop indices
INTEGER :: i01 ! 0 or 1
REAL(wp) :: z1e3w
REAL(wp), DIMENSION(A2D(0),jpk) :: ztrds ! 3D workspace
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwt, zws, ztrdt ! 3D workspace
!!----------------------------------------------------------------------
!
IF( .NOT. ALLOCATED( trdtx ) ) THEN ! allocate trdtra arrays
IF( trd_tra_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trd_tra : unable to allocate arrays' )
avt_evd(:,:,:) = 0._wp
avt_evd(A2D(0),:) = 0._wp
ENDIF
!
i01 = COUNT( (/ PRESENT(pu) .OR. ( ktrd /= jptra_xad .AND. ktrd /= jptra_yad .AND. ktrd /= jptra_zad ) /) )
......@@ -103,13 +104,13 @@ CONTAINS
CASE( jptra_yad ) ; CALL trd_tra_adv( ptrd, pu, ptra, 'Y', trdty, Kmm )
CASE( jptra_zad ) ; CALL trd_tra_adv( ptrd, pu, ptra, 'Z', trdt, Kmm )
CASE( jptra_bbc, & ! qsr, bbc: on temperature only, send to trd_tra_mng
& jptra_qsr ) ; trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
ztrds(:,:,:) = 0._wp
& jptra_qsr ) ; trdt(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
ztrds(A2D(0),:) = 0._wp
CALL trd_tra_mng( trdt, ztrds, ktrd, kt, Kmm )
!!gm Gurvan, verify the jptra_evd trend please !
CASE( jptra_evd ) ; avt_evd(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
CASE( jptra_evd ) ; avt_evd(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
CASE DEFAULT ! other trends: masked trends
trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:) ! mask & store
trdt(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:) ! mask & store
END SELECT
!
ENDIF
......@@ -127,44 +128,44 @@ CONTAINS
CALL trd_tra_mng( trdt , ztrds, ktrd, kt, Kmm )
CASE( jptra_zdfp ) ! diagnose the "PURE" Kz trend (here: just before the swap)
! ! iso-neutral diffusion case otherwise jptra_zdf is "PURE"
ALLOCATE( zwt(jpi,jpj,jpk), zws(jpi,jpj,jpk), ztrdt(jpi,jpj,jpk) )
ALLOCATE( zwt(A2D(0),jpk), zws(A2D(0),jpk), ztrdt(A2D(0),jpk) )
!
zwt(:,:, 1 ) = 0._wp ; zws(:,:, 1 ) = 0._wp ! vertical diffusive fluxes
zwt(:,:,jpk) = 0._wp ; zws(:,:,jpk) = 0._wp
DO jk = 2, jpk
zwt(:,:,jk) = avt(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) ) &
& / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
zws(:,:,jk) = avs(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) ) &
& / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
END DO
zwt(A2D(0), 1 ) = 0._wp ; zws(A2D(0), 1 ) = 0._wp ! vertical diffusive fluxes
zwt(A2D(0),jpk) = 0._wp ; zws(A2D(0),jpk) = 0._wp
DO_3D( 0, 0, 0, 0, 2, jpk )
z1e3w = 1._wp / ( e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) )
zwt(ji,jj,jk) = avt(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_tem,Krhs) - ts(ji,jj,jk,jp_tem,Krhs) ) * z1e3w
zws(ji,jj,jk) = avs(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_sal,Krhs) - ts(ji,jj,jk,jp_sal,Krhs) ) * z1e3w
END_3D
!
ztrdt(:,:,jpk) = 0._wp ; ztrds(:,:,jpk) = 0._wp
DO jk = 1, jpkm1
ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
END DO
ztrdt(A2D(0),jpk) = 0._wp ; ztrds(A2D(0),jpk) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z1e3w = 1._wp / e3w(ji,jj,jk,Kmm)
ztrdt(ji,jj,jk) = ( zwt(ji,jj,jk) - zwt(ji,jj,jk+1) ) * z1e3w
ztrds(ji,jj,jk) = ( zws(ji,jj,jk) - zws(ji,jj,jk+1) ) * z1e3w
END_3D
CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt, Kmm )
!
! ! Also calculate EVD trend at this point.
zwt(:,:,:) = 0._wp ; zws(:,:,:) = 0._wp ! vertical diffusive fluxes
DO jk = 2, jpk
zwt(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) ) &
& / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
zws(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) ) &
& / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
END DO
zwt(A2D(0), :) = 0._wp ; zws(A2D(0), :) = 0._wp ! vertical diffusive fluxes
DO_3D( 0, 0, 0, 0, 2, jpk )
z1e3w = 1._wp / ( e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) )
zwt(ji,jj,jk) = avt_evd(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_tem,Krhs) - ts(ji,jj,jk,jp_tem,Krhs) ) * z1e3w
zws(ji,jj,jk) = avt_evd(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_sal,Krhs) - ts(ji,jj,jk,jp_sal,Krhs) ) * z1e3w
END_3D
!
ztrdt(:,:,jpk) = 0._wp ; ztrds(:,:,jpk) = 0._wp
DO jk = 1, jpkm1
ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
END DO
ztrdt(A2D(0),jpk) = 0._wp ; ztrds(A2D(0),jpk) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z1e3w = 1._wp / e3w(ji,jj,jk,Kmm)
ztrdt(ji,jj,jk) = ( zwt(ji,jj,jk) - zwt(ji,jj,jk+1) ) * z1e3w
ztrds(ji,jj,jk) = ( zws(ji,jj,jk) - zws(ji,jj,jk+1) ) * z1e3w
END_3D
CALL trd_tra_mng( ztrdt, ztrds, jptra_evd, kt, Kmm )
!
DEALLOCATE( zwt, zws, ztrdt )
!
CASE DEFAULT ! other trends: mask and send T & S trends to trd_tra_mng
ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
ztrds(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
CALL trd_tra_mng( trdt, ztrds, ktrd, kt, Kmm )
END SELECT
ENDIF
......@@ -177,7 +178,7 @@ CONTAINS
CASE( jptra_yad ) ; CALL trd_tra_adv( ptrd , pu , ptra, 'Y', ztrds, Kmm )
CASE( jptra_zad ) ; CALL trd_tra_adv( ptrd , pu , ptra, 'Z', ztrds, Kmm )
CASE DEFAULT ! other trends: just masked
ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
ztrds(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
END SELECT
! ! send trend to trd_trc
CALL trd_trc( ztrds, ktra, ktrd, kt, Kmm )
......@@ -199,12 +200,12 @@ CONTAINS
!! k-advective trends = -un. di+1[T] = -( dk+1[fi] - tn dk+1[un] )
!! where fi is the incoming advective flux.
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pf ! advective flux in one direction
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu ! now velocity in one direction
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt ! now or before tracer
CHARACTER(len=1) , INTENT(in ) :: cdir ! X/Y/Z direction
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: ptrd ! advective trend in one direction
INTEGER, INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pf ! advective flux in one direction
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pu ! now velocity in one direction
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! now or before tracer
CHARACTER(len=1) , INTENT(in ) :: cdir ! X/Y/Z direction
REAL(wp), DIMENSION(:,:,:), INTENT( out) :: ptrd ! advective trend in one direction
INTEGER, INTENT(in) :: Kmm ! time level index
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ii, ij, ik ! index shift as function of the direction
......@@ -217,9 +218,7 @@ CONTAINS
END SELECT
!
! ! set to zero uncomputed values
ptrd(jpi,:,:) = 0._wp ; ptrd(1,:,:) = 0._wp
ptrd(:,jpj,:) = 0._wp ; ptrd(:,1,:) = 0._wp
ptrd(:,:,jpk) = 0._wp
ptrd(:,:,:) = 0._wp
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! advective trend
ptrd(ji,jj,jk) = - ( pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik) &
......@@ -248,7 +247,7 @@ CONTAINS
! ! 3D output of tracers trends using IOM interface
IF( ln_tra_trd ) CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt, Kmm )
! ! Integral Constraints Properties for tracers trends !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! ! Integral Constraints Properties for tracers trends
IF( ln_glo_trd ) CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt, Kmm )
! ! Potential ENergy trends
......@@ -305,8 +304,8 @@ CONTAINS
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! time level index
!!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ikbu, ikbv ! local integers
INTEGER :: ji, jj
REAL(wp) :: z1e3
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2dx, z2dy ! 2D workspace
!!----------------------------------------------------------------------
!
......@@ -329,9 +328,12 @@ CONTAINS
CASE( jptra_zad ) ; CALL iom_put( "ttrd_zad" , ptrdx ) ! z- vertical advection
CALL iom_put( "strd_zad" , ptrdy )
IF( ln_linssh ) THEN ! cst volume : adv flux through z=0 surface
ALLOCATE( z2dx(jpi,jpj), z2dy(jpi,jpj) )
z2dx(:,:) = ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) / e3t(:,:,1,Kmm)
z2dy(:,:) = ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) / e3t(:,:,1,Kmm)
ALLOCATE( z2dx(A2D(0)), z2dy(A2D(0)) )
DO_2D( 0, 0, 0, 0 )
z1e3 = ww(ji,jj,1) / e3t(ji,jj,1,Kmm)
z2dx(ji,jj) = ts(ji,jj,1,jp_tem,Kmm) * z1e3
z2dy(ji,jj) = ts(ji,jj,1,jp_sal,Kmm) * z1e3
END_2D
CALL iom_put( "ttrd_sad", z2dx )
CALL iom_put( "strd_sad", z2dy )
DEALLOCATE( z2dx, z2dy )
......
src/OCE/TRD/trdvor.F90
View file @ ee6ce8c5
......@@ -68,9 +68,9 @@ CONTAINS
!!----------------------------------------------------------------------------
!! *** ROUTINE trd_vor_alloc ***
!!----------------------------------------------------------------------------
ALLOCATE( vor_avr (jpi,jpj) , vor_avrb(jpi,jpj) , vor_avrbb (jpi,jpj) , &
& vor_avrbn (jpi,jpj) , rotot (jpi,jpj) , vor_avrtot(jpi,jpj) , &
& vor_avrres(jpi,jpj) , vortrd (jpi,jpj,jpltot_vor) , &
ALLOCATE( vor_avr (A2D(0)) , vor_avrb(A2D(0)) , vor_avrbb (A2D(0)) , &
& vor_avrbn (A2D(0)) , rotot (A2D(0)) , vor_avrtot(A2D(0)) , &
& vor_avrres(A2D(0)) , vortrd (A2D(0),jpltot_vor) , &
& ndexvor1 (jpi*jpj) , STAT= trd_vor_alloc )
!
CALL mpp_sum ( 'trdvor', trd_vor_alloc )
......@@ -105,9 +105,9 @@ CONTAINS
CASE( jpdyn_zad ) ; CALL trd_vor_zint( putrd, pvtrd, jpvor_zad, Kmm ) ! Vertical Advection
CASE( jpdyn_spg ) ; CALL trd_vor_zint( putrd, pvtrd, jpvor_spg, Kmm ) ! Surface Pressure Grad.
CASE( jpdyn_zdf ) ! Vertical Diffusion
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! wind stress trends
ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
DO_2D( 0, 1, 0, 1 ) ! wind stress trends
ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utauU(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtauV(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
END_2D
CALL trd_vor_zint( putrd, pvtrd, jpvor_zdf, Kmm ) ! zdf trend including surf./bot. stresses
CALL trd_vor_zint( ztswu, ztswv, jpvor_swf, Kmm ) ! surface wind stress
......@@ -165,7 +165,7 @@ CONTAINS
SELECT CASE( ktrd )
!
CASE( jpvor_bfr ) ! bottom friction
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 1, 0, 1 )
ikbu = mbkv(ji,jj)
ikbv = mbkv(ji,jj)
zudpvor(ji,jj) = putrdvor(ji,jj) * e3u(ji,jj,ikbu,Kmm) * e1u(ji,jj) * umask(ji,jj,ikbu)
......@@ -173,14 +173,18 @@ CONTAINS
END_2D
!
CASE( jpvor_swf ) ! wind stress
zudpvor(:,:) = putrdvor(:,:) * e3u(:,:,1,Kmm) * e1u(:,:) * umask(:,:,1)
zvdpvor(:,:) = pvtrdvor(:,:) * e3v(:,:,1,Kmm) * e2v(:,:) * vmask(:,:,1)
DO_2D( 0, 1, 0, 1 )
zudpvor(ji,jj) = putrdvor(ji,jj) * e3u(ji,jj,1,Kmm) * e1u(ji,jj) * umask(ji,jj,1)
zvdpvor(ji,jj) = pvtrdvor(ji,jj) * e3v(ji,jj,1,Kmm) * e2v(ji,jj) * vmask(ji,jj,1)
END_2D
!
END SELECT
! Average except for Beta.V
zudpvor(:,:) = zudpvor(:,:) * r1_hu(:,:,Kmm)
zvdpvor(:,:) = zvdpvor(:,:) * r1_hv(:,:,Kmm)
DO_2D( 0, 1, 0, 1 )
zudpvor(ji,jj) = zudpvor(ji,jj) * r1_hu(ji,jj,Kmm)
zvdpvor(ji,jj) = zvdpvor(ji,jj) * r1_hv(ji,jj,Kmm)
END_2D
! Curl
DO_2D( 0, 0, 0, 0 )
......@@ -238,10 +242,11 @@ CONTAINS
! I vertical integration of 3D trends
! =====================================
! putrdvor and pvtrdvor terms
DO jk = 1,jpk
zudpvor(:,:) = zudpvor(:,:) + putrdvor(:,:,jk) * e3u(:,:,jk,Kmm) * e1u(:,:) * umask(:,:,jk)
zvdpvor(:,:) = zvdpvor(:,:) + pvtrdvor(:,:,jk) * e3v(:,:,jk,Kmm) * e2v(:,:) * vmask(:,:,jk)
END DO
zudpvor(:,:) = 0._wp ; zvdpvor(:,:) = 0._wp
DO_3D( 0, 1, 0, 1, 1, jpk )
zudpvor(ji,jj) = zudpvor(ji,jj) + putrdvor(ji,jj,jk) * e3u(ji,jj,jk,Kmm) * e1u(ji,jj) * umask(ji,jj,jk)
zvdpvor(ji,jj) = zvdpvor(ji,jj) + pvtrdvor(ji,jj,jk) * e3v(ji,jj,jk,Kmm) * e2v(ji,jj) * vmask(ji,jj,jk)
END_3D
! Planetary vorticity: 2nd computation (Beta.V term) store the vertical sum
! as Beta.V term need intergration, not average
......@@ -254,8 +259,10 @@ CONTAINS
ENDIF
!
! Average
zudpvor(:,:) = zudpvor(:,:) * r1_hu(:,:,Kmm)
zvdpvor(:,:) = zvdpvor(:,:) * r1_hv(:,:,Kmm)
DO_2D( 0, 1, 0, 1 )
zudpvor(ji,jj) = zudpvor(ji,jj) * r1_hu(ji,jj,Kmm)
zvdpvor(ji,jj) = zvdpvor(ji,jj) * r1_hv(ji,jj,Kmm)
END_2D
!
! Curl
DO_2D( 0, 0, 0, 0 )
......@@ -368,10 +375,6 @@ CONTAINS
zmean = 1._wp / REAL( nmoydpvor, wp )
vor_avrres(:,:) = vor_avrtot(:,:) - rotot(:,:) / zmean
! Boundary conditions
CALL lbc_lnk( 'trdvor', vor_avrtot, 'F', 1.0_wp , vor_avrres, 'F', 1.0_wp )
! III.3 time evolution array swap
! ------------------------------
vor_avrbb(:,:) = vor_avrb(:,:)
......