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
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 431 additions and 318 deletions
src/ICE/iceupdate.F90
View file @ 32d63bb7
......@@ -314,7 +314,7 @@ CONTAINS
!!
!! ** Action : * at each ice time step (every nn_fsbc time step):
!! - compute the modulus of ice-ocean relative velocity
!! (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
!! (*rho*Cd) at T-point (C-grid)
!! tmod_io = rhoco * | U_ice-U_oce |
!! - update the modulus of stress at ocean surface
!! taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce |
......@@ -325,19 +325,19 @@ CONTAINS
!!
!! NB: - ice-ocean rotation angle no more allowed
!! - here we make an approximation: taum is only computed every ice time step
!! This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
!! to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
!! This avoids mutiple average to pass from U,V grids to T grids
!! taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
!!
!! ** Outputs : - utau, vtau : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
!! ** Outputs : - utau, vtau : surface ocean i- and j-stress (T-pts) updated with ice-ocean fluxes
!! - taum : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
!!---------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! ocean time-step index
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zat_u, zutau_ice, zu_t, zmodt ! local scalar
REAL(wp) :: zat_v, zvtau_ice, zv_t, zrhoco ! - -
REAL(wp) :: zflagi ! - -
REAL(wp) :: zutau_ice, zu_t, zmodt ! local scalar
REAL(wp) :: zvtau_ice, zv_t, zrhoco ! - -
REAL(wp) :: zflagi ! - -
!!---------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('iceupdate')
......@@ -352,8 +352,8 @@ CONTAINS
IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
DO_2D( 0, 0, 0, 0 ) !* update the modulus of stress at ocean surface (T-point)
! ! 2*(U_ice-U_oce) at T-point
zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)
zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)
zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj) ! u_oce = ssu_m
zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) ! v_oce = ssv_m
! ! |U_ice-U_oce|^2
zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t )
! ! update the ocean stress modulus
......@@ -377,19 +377,14 @@ CONTAINS
ENDIF
!
DO_2D( 0, 0, 0, 0 ) !* update the stress WITHOUT an ice-ocean rotation angle
! ice area at u and v-points
zat_u = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji+1,jj ) * tmask(ji+1,jj ,1) ) &
& / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji+1,jj ,1) )
zat_v = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji ,jj+1 ) * tmask(ji ,jj+1,1) ) &
& / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji ,jj+1,1) )
! ! linearized quadratic drag formulation
zutau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
zvtau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
zutau_ice = 0.5_wp * tmod_io(ji,jj) * ( u_ice(ji,jj) + u_ice(ji-1,jj) - pu_oce(ji,jj) - pu_oce(ji-1,jj) )
zvtau_ice = 0.5_wp * tmod_io(ji,jj) * ( v_ice(ji,jj) + v_ice(ji,jj-1) - pv_oce(ji,jj) - pv_oce(ji,jj-1) )
! ! stresses at the ocean surface
utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
utau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * utau_oce(ji,jj) + at_i(ji,jj) * zutau_ice
vtau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * vtau_oce(ji,jj) + at_i(ji,jj) * zvtau_ice
END_2D
CALL lbc_lnk( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp ) ! lateral boundary condition
CALL lbc_lnk( 'iceupdate', utau, 'T', -1.0_wp, vtau, 'T', -1.0_wp ) ! lateral boundary condition
!
IF( ln_timing ) CALL timing_stop('iceupdate')
!
......
src/OCE/CRS/crsfld.F90
View file @ 32d63bb7
......@@ -211,8 +211,8 @@ CONTAINS
! sbc fields
CALL crs_dom_ope( ssh(:,:,Kmm) , 'VOL', 'T', tmask, sshn_crs , p_e12=e1e2t, p_e3=ze3t , psgn=1.0_wp )
CALL crs_dom_ope( utau , 'SUM', 'U', umask, utau_crs , p_e12=e2u , p_surf_crs=e2u_crs , psgn=1.0_wp )
CALL crs_dom_ope( vtau , 'SUM', 'V', vmask, vtau_crs , p_e12=e1v , p_surf_crs=e1v_crs , psgn=1.0_wp )
CALL crs_dom_ope( utau , 'SUM', 'T', tmask, utau_crs , p_e12=e2t , p_surf_crs=e2t_crs , psgn=1.0_wp ) !clem tau: check psgn ??
CALL crs_dom_ope( vtau , 'SUM', 'T', tmask, vtau_crs , p_e12=e1t , p_surf_crs=e1t_crs , psgn=1.0_wp ) !
CALL crs_dom_ope( wndm , 'SUM', 'T', tmask, wndm_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0_wp )
CALL crs_dom_ope( rnf , 'MAX', 'T', tmask, rnf_crs , psgn=1.0_wp )
CALL crs_dom_ope( qsr , 'SUM', 'T', tmask, qsr_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0_wp )
......
src/OCE/DIA/diadetide.F90
View file @ 32d63bb7
......@@ -5,11 +5,11 @@ MODULE diadetide
!!======================================================================
!! History : ! 2019 (S. Mueller)
!!----------------------------------------------------------------------
USE par_oce , ONLY : wp, jpi, jpj
USE in_out_manager , ONLY : lwp, numout
USE iom , ONLY : iom_put
USE dom_oce , ONLY : rn_Dt, nsec_day
USE phycst , ONLY : rpi
USE par_oce
USE in_out_manager
USE iom
USE dom_oce
USE phycst
USE tide_mod
#if defined key_xios
USE xios
......@@ -24,6 +24,8 @@ MODULE diadetide
PUBLIC :: dia_detide_init, dia_detide
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2019)
!! $Id$
......@@ -90,9 +92,9 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt
REAL(wp), DIMENSION(jpi,jpj) :: zwght_2D
REAL(wp), DIMENSION(A2D(0)) :: zwght_2D
REAL(wp) :: zwght, ztmp
INTEGER :: jn
INTEGER :: ji, jj, jn
! Compute detiding weight at the current time-step; the daily total weight
! is one, and the daily summation of a diagnosed field multiplied by this
......@@ -104,7 +106,10 @@ CONTAINS
zwght = zwght + 1.0_wp / REAL( ndiadetide, KIND=wp )
END IF
END DO
zwght_2D(:,:) = zwght
DO_2D( 0, 0, 0, 0 )
zwght_2D(ji,jj) = zwght
END_2D
CALL iom_put( "diadetide_weight", zwght_2D)
END SUBROUTINE dia_detide
......
src/OCE/DIA/diahsb.F90
View file @ 32d63bb7
......@@ -50,6 +50,7 @@ MODULE diahsb
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmask_ini
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -82,30 +83,61 @@ CONTAINS
REAL(wp) :: z_frc_trd_v ! - -
REAL(wp) :: z_wn_trd_t , z_wn_trd_s ! - -
REAL(wp) :: z_ssh_hc , z_ssh_sc ! - -
REAL(wp), DIMENSION(jpi,jpj,13) :: ztmp
REAL(wp), DIMENSION(jpi,jpj,jpkm1,4) :: ztmpk
REAL(wp), DIMENSION(A2D(0),13) :: ztmp
REAL(wp), DIMENSION(A2D(0),jpkm1,4) :: ztmpk
REAL(wp), DIMENSION(17) :: zbg
!!---------------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('dia_hsb')
!
ztmp (:,:,:) = 0._wp ! should be better coded
ztmpk(:,:,:,:) = 0._wp ! should be better coded
!
ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ;
ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ;
DO_2D( 0, 0, 0, 0 )
ztmp (ji,jj,:) = 0._wp ! should be better coded
ztmpk(ji,jj,:,:) = 0._wp ! should be better coded
!
ts(ji,jj,:,1,Kmm) = ts(ji,jj,:,1,Kmm) * tmask(ji,jj,:)
ts(ji,jj,:,1,Kbb) = ts(ji,jj,:,1,Kbb) * tmask(ji,jj,:)
!
ts(ji,jj,:,2,Kmm) = ts(ji,jj,:,2,Kmm) * tmask(ji,jj,:)
ts(ji,jj,:,2,Kbb) = ts(ji,jj,:,2,Kbb) * tmask(ji,jj,:)
END_2D
!
! ------------------------- !
! 1 - Trends due to forcing !
! ------------------------- !
! prepare trends
ztmp(:,:,1) = - r1_rho0 * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) ) * surf(:,:) ! volume
ztmp(:,:,2) = sbc_tsc(:,:,jp_tem) * surf(:,:) ! heat
ztmp(:,:,3) = sbc_tsc(:,:,jp_sal) * surf(:,:) ! salt
IF( ln_rnf ) ztmp(:,:,4) = rnf_tsc(:,:,jp_tem) * surf(:,:) ! runoff temp
IF( ln_rnf_sal ) ztmp(:,:,5) = rnf_tsc(:,:,jp_sal) * surf(:,:) ! runoff salt
IF( ln_isf ) ztmp(:,:,6) = ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) ! isf temp
IF( ln_traqsr ) ztmp(:,:,7) = r1_rho0_rcp * qsr(:,:) * surf(:,:) ! penetrative solar radiation
IF( ln_trabbc ) ztmp(:,:,8) = qgh_trd0(:,:) * surf(:,:) ! geothermal heat
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,1) = - r1_rho0 * ( emp(ji,jj) & ! volume
& - rnf(ji,jj) &
& - fwfisf_cav(ji,jj) &
& - fwfisf_par(ji,jj) ) * surf(ji,jj)
ztmp(ji,jj,2) = sbc_tsc(ji,jj,jp_tem) * surf(ji,jj) ! heat
ztmp(ji,jj,3) = sbc_tsc(ji,jj,jp_sal) * surf(ji,jj) ! salt
END_2D
IF( ln_rnf ) THEN
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,4) = rnf_tsc(ji,jj,jp_tem) * surf(ji,jj) ! runoff temp
END_2D
END IF
IF( ln_rnf_sal ) THEN
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,5) = rnf_tsc(ji,jj,jp_sal) * surf(ji,jj) ! runoff salt
END_2D
END IF
IF( ln_isf ) THEN
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,6) = ( risf_cav_tsc(ji,jj,jp_tem) &
& + risf_par_tsc(ji,jj,jp_tem) ) * surf(ji,jj) ! isf temp
END_2D
END IF
IF( ln_traqsr ) THEN
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,7) = r1_rho0_rcp * qsr(ji,jj) * surf(ji,jj) ! penetrative solar radiation
END_2D
END IF
IF( ln_trabbc ) THEN
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,8) = qgh_trd0(ji,jj) * surf(ji,jj) ! geothermal heat
END_2D
END IF
!
IF( ln_linssh ) THEN ! Advection flux through fixed surface (z=0)
IF( ln_isfcav ) THEN
......@@ -116,8 +148,10 @@ CONTAINS
END DO
END DO
ELSE
ztmp(:,:,9 ) = - surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
ztmp(:,:,10) = - surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,9 ) = - surf(ji,jj) * ww(ji,jj,1) * ts(ji,jj,1,jp_tem,Kbb)
ztmp(ji,jj,10) = - surf(ji,jj) * ww(ji,jj,1) * ts(ji,jj,1,jp_sal,Kbb)
END_2D
END IF
ENDIF
......@@ -152,7 +186,9 @@ CONTAINS
! glob_sum is needed because you keep only the interior domain to compute the sum (iscpl)
!
! ! volume variation (calculated with ssh)
ztmp(:,:,11) = surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:)
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,11) = surf(ji,jj)*ssh(ji,jj,Kmm) - surf_ini(ji,jj)*ssh_ini(ji,jj)
END_2D
! ! heat & salt content variation (associated with ssh)
IF( ln_linssh ) THEN ! linear free surface case
......@@ -164,8 +200,10 @@ CONTAINS
END DO
END DO
ELSE ! no under ice-shelf seas
ztmp(:,:,12) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) )
ztmp(:,:,13) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) )
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,12) = surf(ji,jj) * ( ts(ji,jj,1,jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )
ztmp(ji,jj,13) = surf(ji,jj) * ( ts(ji,jj,1,jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )
END_2D
END IF
ENDIF
......@@ -185,19 +223,27 @@ CONTAINS
! glob_sum is needed because you keep only the interior domain to compute the sum (iscpl)
!
DO jk = 1, jpkm1 ! volume
ztmpk(:,:,jk,1) = surf (:,:) * e3t(:,:,jk,Kmm)*tmask(:,:,jk) &
& - surf_ini(:,:) * e3t_ini(:,:,jk )*tmask_ini(:,:,jk)
DO_2D( 0, 0, 0, 0 )
ztmpk(ji,jj,jk,1) = surf (ji,jj) * e3t(ji,jj,jk,Kmm)*tmask(ji,jj,jk) &
& - surf_ini(ji,jj) * e3t_ini(ji,jj,jk )*tmask_ini(ji,jj,jk)
END_2D
END DO
DO jk = 1, jpkm1 ! heat
ztmpk(:,:,jk,2) = ( surf (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm) &
& - surf_ini(:,:) * hc_loc_ini(:,:,jk) )
DO_2D( 0, 0, 0, 0 )
ztmpk(ji,jj,jk,2) = ( surf (ji,jj) * e3t(ji,jj,jk,Kmm)*ts(ji,jj,jk,jp_tem,Kmm) &
& - surf_ini(ji,jj) * hc_loc_ini(ji,jj,jk) )
END_2D
END DO
DO jk = 1, jpkm1 ! salt
ztmpk(:,:,jk,3) = ( surf (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm) &
& - surf_ini(:,:) * sc_loc_ini(:,:,jk) )
DO_2D( 0, 0, 0, 0 )
ztmpk(ji,jj,jk,3) = ( surf (ji,jj) * e3t(ji,jj,jk,Kmm)*ts(ji,jj,jk,jp_sal,Kmm) &
& - surf_ini(ji,jj) * sc_loc_ini(ji,jj,jk) )
END_2D
END DO
DO jk = 1, jpkm1 ! total ocean volume
ztmpk(:,:,jk,4) = surf(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
DO_2D( 0, 0, 0, 0 )
ztmpk(ji,jj,jk,4) = surf(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
END_2D
END DO
! global sum
......@@ -315,14 +361,18 @@ CONTAINS
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' dia_hsb_rst : initialise hsb at initial state '
IF(lwp) WRITE(numout,*)
surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:) ! initial ocean surface
ssh_ini(:,:) = ssh(:,:,Kmm) ! initial ssh
DO_2D( 0, 0, 0, 0 )
surf_ini(ji,jj) = e1e2t(ji,jj) * tmask_i(ji,jj) ! initial ocean surface
ssh_ini(ji,jj) = ssh(ji,jj,Kmm) ! initial ssh
END_2D
DO jk = 1, jpk
! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
e3t_ini (:,:,jk) = e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial vertical scale factors
tmask_ini (:,:,jk) = tmask(:,:,jk) ! initial mask
hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial heat content
sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial salt content
DO_2D( 0, 0, 0, 0 )
e3t_ini (ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) ! initial vertical scale factors
tmask_ini (ji,jj,jk) = tmask(ji,jj,jk) ! initial mask
hc_loc_ini(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) ! initial heat content
sc_loc_ini(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) ! initial salt content
END_2D
END DO
frc_v = 0._wp ! volume trend due to forcing
frc_t = 0._wp ! heat content - - - -
......@@ -334,13 +384,15 @@ CONTAINS
ssh_hc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) ! initial heat content in ssh
ssh_sc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) ! initial salt content in ssh
END DO
END DO
ELSE
ssh_hc_loc_ini(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) ! initial heat content in ssh
ssh_sc_loc_ini(:,:) = ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) ! initial salt content in ssh
END DO
ELSE
DO_2D( 0, 0, 0, 0 )
ssh_hc_loc_ini(ji,jj) = ts(ji,jj,1,jp_tem,Kmm) * ssh(ji,jj,Kmm) ! initial heat content in ssh
ssh_sc_loc_ini(ji,jj) = ts(ji,jj,1,jp_sal,Kmm) * ssh(ji,jj,Kmm) ! initial salt content in ssh
END_2D
END IF
frc_wn_t = 0._wp ! initial heat content misfit due to free surface
frc_wn_s = 0._wp ! initial salt content misfit due to free surface
frc_wn_t = 0._wp ! initial heat content misfit due to free surface
frc_wn_s = 0._wp ! initial salt content misfit due to free surface
ENDIF
ENDIF
!
......@@ -388,6 +440,7 @@ CONTAINS
INTEGER, INTENT(in) :: Kmm ! time level index
!
INTEGER :: ierror, ios ! local integer
INTEGER :: ji, jj ! loop index
!!
NAMELIST/namhsb/ ln_diahsb
!!----------------------------------------------------------------------
......@@ -427,7 +480,10 @@ CONTAINS
! ----------------------------------------------- !
! 2 - Time independant variables and file opening !
! ----------------------------------------------- !
surf(:,:) = e1e2t(:,:) * tmask_i(:,:) ! masked surface grid cell area
DO_2D( 0, 0, 0, 0 )
surf(ji,jj) = e1e2t(ji,jj) * tmask_i(ji,jj) ! masked surface grid cell area
END_2D
surf_tot = glob_sum( 'diahsb', surf(:,:) ) ! total ocean surface area
IF( ln_bdy ) CALL ctl_warn( 'dia_hsb_init: heat/salt budget does not consider open boundary fluxes' )
......
src/OCE/DIA/diahth.F90
View file @ 32d63bb7
......@@ -86,22 +86,22 @@ CONTAINS
INTEGER, INTENT( in ) :: kt ! ocean time-step index
INTEGER, INTENT( in ) :: Kmm ! ocean time level index
!!
INTEGER :: ji, jj, jk ! dummy loop arguments
REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth
REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth
REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth
REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop
REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace
REAL(wp), DIMENSION(jpi,jpj) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2
REAL(wp), DIMENSION(jpi,jpj) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2
REAL(wp), DIMENSION(jpi,jpj) :: zrho10_3 ! MLD: rho = rho10m + zrho3
REAL(wp), DIMENSION(jpi,jpj) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
REAL(wp), DIMENSION(jpi,jpj) :: ztinv ! max of temperature inversion
REAL(wp), DIMENSION(jpi,jpj) :: zdepinv ! depth of temperature inversion
REAL(wp), DIMENSION(jpi,jpj) :: zrho0_3 ! MLD rho = rho(surf) = 0.03
REAL(wp), DIMENSION(jpi,jpj) :: zrho0_1 ! MLD rho = rho(surf) = 0.01
REAL(wp), DIMENSION(jpi,jpj) :: zmaxdzT ! max of dT/dz
REAL(wp), DIMENSION(jpi,jpj) :: zdelr ! delta rho equivalent to deltaT = 0.2
INTEGER :: ji, jj, jk ! dummy loop arguments
REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth
REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth
REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth
REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop
REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace
REAL(wp), DIMENSION(A2D(0)) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2
REAL(wp), DIMENSION(A2D(0)) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2
REAL(wp), DIMENSION(A2D(0)) :: zrho10_3 ! MLD: rho = rho10m + zrho3
REAL(wp), DIMENSION(A2D(0)) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
REAL(wp), DIMENSION(A2D(0)) :: ztinv ! max of temperature inversion
REAL(wp), DIMENSION(A2D(0)) :: zdepinv ! depth of temperature inversion
REAL(wp), DIMENSION(A2D(0)) :: zrho0_3 ! MLD rho = rho(surf) = 0.03
REAL(wp), DIMENSION(A2D(0)) :: zrho0_1 ! MLD rho = rho(surf) = 0.01
REAL(wp), DIMENSION(A2D(0)) :: zmaxdzT ! max of dT/dz
REAL(wp), DIMENSION(A2D(0)) :: zdelr ! delta rho equivalent to deltaT = 0.2
!!----------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('dia_hth')
......@@ -131,7 +131,7 @@ CONTAINS
IF( iom_use( 'mlddzt' ) ) zmaxdzT(:,:) = 0._wp
IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) &
& .OR. iom_use( 'mldr10_3' ) .OR. iom_use( 'pycndep' ) ) THEN
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)
hth (ji,jj) = zztmp
zabs2 (ji,jj) = zztmp
......@@ -142,7 +142,7 @@ CONTAINS
ENDIF
IF( iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
IF( nla10 > 1 ) THEN
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)
zrho0_3(ji,jj) = zztmp
zrho0_1(ji,jj) = zztmp
......@@ -157,7 +157,7 @@ CONTAINS
! MLD: rho = rho(1) + zrho3 !
! MLD: rho = rho(1) + zrho1 !
! ------------------------------------------------------------- !
DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) ! loop from bottom to 2
DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! loop from bottom to 2
!
zzdep = gdepw(ji,jj,jk,Kmm)
zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) &
......@@ -189,7 +189,7 @@ CONTAINS
!
! Preliminary computation
! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
IF( tmask(ji,jj,nla10) == 1. ) THEN
zu = 1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80 * ts(ji,jj,nla10,jp_sal,Kmm) &
& - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) &
......@@ -213,7 +213,7 @@ CONTAINS
! temperature inversion: max( 0, max of tn - tn(10m) ) !
! depth of temperature inversion !
! ------------------------------------------------------------- !
DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 ) ! loop from bottom to nlb10
DO_3DS( 0, 0, 0, 0, jpkm1, nlb10, -1 ) ! loop from bottom to nlb10
!
zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1)
!
......@@ -305,13 +305,16 @@ CONTAINS
!
INTEGER :: ji, jj, jk, iid
REAL(wp) :: zztmp, zzdep
INTEGER, DIMENSION(jpi,jpj) :: iktem
INTEGER, DIMENSION(A2D(0)) :: iktem
! --------------------------------------- !
! search deepest level above ptem !
! --------------------------------------- !
iktem(:,:) = 1
DO_3D( 1, 1, 1, 1, 1, jpkm1 ) ! beware temperature is not always decreasing with depth => loop from top to bottom
DO_2D( 0, 0, 0, 0 )
iktem(ji,jj) = 1
END_2D
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! beware temperature is not always decreasing with depth => loop from top to bottom
zztmp = ts(ji,jj,jk,jp_tem,Kmm)
IF( zztmp >= ptem ) iktem(ji,jj) = jk
END_3D
......@@ -319,7 +322,7 @@ CONTAINS
! ------------------------------- !
! Depth of ptem isotherm !
! ------------------------------- !
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
!
zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! depth of the ocean bottom
!
......@@ -346,18 +349,29 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc
!
INTEGER :: ji, jj, jk, ik
REAL(wp), DIMENSION(jpi,jpj) :: zthick
INTEGER , DIMENSION(jpi,jpj) :: ilevel
REAL(wp), DIMENSION(A2D(0)) :: zthick
INTEGER , DIMENSION(A2D(0)) :: ilevel
! surface boundary condition
IF( .NOT. ln_linssh ) THEN ; zthick(:,:) = 0._wp ; phtc(:,:) = 0._wp
ELSE ; zthick(:,:) = ssh(:,:,Kmm) ; phtc(:,:) = pt(:,:,1) * ssh(:,:,Kmm) * tmask(:,:,1)
IF( .NOT. ln_linssh ) THEN
DO_2D( 0, 0, 0, 0 )
zthick(ji,jj) = 0._wp
phtc (ji,jj) = 0._wp
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
zthick(ji,jj) = ssh(ji,jj,Kmm)
phtc (ji,jj) = pt(ji,jj,1) * ssh(ji,jj,Kmm) * tmask(ji,jj,1)
END_2D
ENDIF
!
ilevel(:,:) = 1
DO_3D( 1, 1, 1, 1, 1, jpkm1 )
DO_2D( 0, 0, 0, 0 )
ilevel(ji,jj) = 1
END_2D
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
IF( ( gdepw(ji,jj,jk+1,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
ilevel(ji,jj) = jk+1
zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm)
......@@ -365,7 +379,7 @@ CONTAINS
ENDIF
END_3D
!
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
ik = ilevel(ji,jj)
IF( tmask(ji,jj,ik) == 1 ) THEN
zthick(ji,jj) = MIN ( gdepw(ji,jj,ik+1,Kmm), pdep ) - zthick(ji,jj) ! remaining thickness to reach dephw pdep
......
src/OCE/DIA/diamlr.F90
View file @ 32d63bb7
......@@ -6,7 +6,7 @@ MODULE diamlr
!! History : 4.0 ! 2019 (S. Mueller) Original code
!!----------------------------------------------------------------------
USE par_oce , ONLY : wp, jpi, jpj
USE par_oce , ONLY : wp, jpi, jpj, ntsi, ntei, ntsj, ntej
USE phycst , ONLY : rpi
USE dom_oce , ONLY : adatrj
USE tide_mod
......@@ -407,8 +407,9 @@ CONTAINS
!! ** Purpose : update time used in multiple-linear-regression analysis
!!
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: zadatrj2d
REAL(wp), DIMENSION(A2D(0)) :: zadatrj2d
!!----------------------------------------------------------------------
INTEGER :: ji, jj
IF( ln_timing ) CALL timing_start('dia_mlr')
......@@ -417,7 +418,9 @@ CONTAINS
!
! A 2-dimensional field of constant value is sent, and subsequently used directly
! or transformed to a scalar or a constant 3-dimensional field as required.
zadatrj2d(:,:) = adatrj*86400.0_wp
DO_2D( 0, 0, 0, 0 )
zadatrj2d(ji,jj) = adatrj*86400.0_wp
END_2D
IF ( iom_use('diamlr_time') ) CALL iom_put('diamlr_time', zadatrj2d)
!
IF( ln_timing ) CALL timing_stop('dia_mlr')
......
src/OCE/DIA/diaptr.F90
View file @ 32d63bb7
This diff is collapsed.
src/OCE/DIA/diawri.F90
View file @ 32d63bb7
......@@ -135,8 +135,8 @@ CONTAINS
ENDIF
! initialize arrays
z2d(:,:) = 0._wp
z3d(:,:,:) = 0._wp
z2d(A2D(0)) = 0._wp
z3d(A2D(0),:) = 0._wp
! Output of initial vertical scale factor
CALL iom_put("e3t_0", e3t_0(:,:,:) )
......@@ -868,7 +868,11 @@ CONTAINS
CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
#endif
CALL histdef( nid_T, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
& jpi, jpj, nh_T, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
CALL histdef( nid_T, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
& jpi, jpj, nh_T, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
!
CALL histend( nid_T, snc4chunks=snc4set )
! !!! nid_U : 3D
......@@ -878,10 +882,7 @@ CONTAINS
CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
ENDIF
! !!! nid_U : 2D
CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
!
CALL histend( nid_U, snc4chunks=snc4set )
! !!! nid_V : 3D
......@@ -891,10 +892,7 @@ CONTAINS
CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
ENDIF
! !!! nid_V : 2D
CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
!
CALL histend( nid_V, snc4chunks=snc4set )
! !!! nid_W : 3D
......@@ -1066,12 +1064,12 @@ CONTAINS
CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
#endif
CALL histwrite( nid_T, "sozotaux", it, utau , ndim_hT, ndex_hT ) ! i-wind stress
CALL histwrite( nid_T, "sometauy", it, vtau , ndim_hT, ndex_hT ) ! j-wind stress
CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm) , ndim_U , ndex_U ) ! i-current
CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm) , ndim_V , ndex_V ) ! j-current
CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
IF( ln_zad_Aimp ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
......
src/OCE/DOM/dom_oce.F90
View file @ 32d63bb7
......@@ -27,6 +27,8 @@ MODULE dom_oce
PUBLIC dom_oce_alloc ! Called from nemogcm.F90
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! time & space domain namelist
!! ----------------------------
......@@ -194,6 +196,7 @@ MODULE dom_oce
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: top first wet T-, U-, V-, F-level (ISF)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: smask0 !: surface mask at T-pts on inner domain
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, wmask, fmask !: land/ocean mask at T-, U-, V-, W- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wumask, wvmask !: land/ocean mask at WU- and WV-pts
......@@ -326,7 +329,7 @@ CONTAINS
ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(ii) )
!
ii = ii+1
ALLOCATE( tmask_i(jpi,jpj) , &
ALLOCATE( tmask_i(jpi,jpj) , smask0(A2D(0)) , &
& ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , &
& mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , mbkf(jpi,jpj) , STAT=ierr(ii) )
!
......
src/OCE/DOM/dommsk.F90
View file @ 32d63bb7
......@@ -144,7 +144,8 @@ CONTAINS
tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
END_3D
ENDIF
smask0(A2D(0)) = tmask(A2D(0),1)
! Ocean/land mask at u-, v-, and f-points (computed from tmask)
! ----------------------------------------
! NB: at this point, fmask is designed for free slip lateral boundary condition
......
src/OCE/DOM/domtile.F90
View file @ 32d63bb7
......@@ -52,16 +52,6 @@ CONTAINS
!!----------------------------------------------------------------------
IF( ln_tile .AND. nn_hls /= 2 ) CALL ctl_stop('dom_tile_init: Tiling is only supported for nn_hls = 2')
ntile = 0 ! Initialise to full domain
nijtile = 1
ntsi = Nis0
ntsj = Njs0
ntei = Nie0
ntej = Nje0
nthl = 0
nthr = 0
nthb = 0
ntht = 0
l_istiled = .FALSE.
IF( ln_tile ) THEN ! Calculate tile domain indices
......
src/OCE/DYN/dynatf.F90
View file @ 32d63bb7
......@@ -331,7 +331,7 @@ CONTAINS
ALLOCATE(zutau(jpi,jpj))
DO_2D( 0, 0, 0, 0 )
jk = miku(ji,jj)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( puu(ji-1,jj,jk,Kaa) + puu(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "utau", zutau(:,:) )
DEALLOCATE(zutau)
......@@ -345,7 +345,7 @@ CONTAINS
ALLOCATE(zvtau(jpi,jpj))
DO_2D( 0, 0, 0, 0 )
jk = mikv(ji,jj)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( pvv(ji,jj-1,jk,Kaa) + pvv(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "vtau", zvtau(:,:) )
DEALLOCATE(zvtau)
......
src/OCE/DYN/dynatf_qco.F90
View file @ 32d63bb7
......@@ -248,7 +248,7 @@ CONTAINS
ALLOCATE(zutau(jpi,jpj))
DO_2D( 0, 0, 0, 0 )
jk = miku(ji,jj)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( puu(ji-1,jj,jk,Kaa) + puu(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "utau", zutau(:,:) )
DEALLOCATE(zutau)
......@@ -262,7 +262,7 @@ CONTAINS
ALLOCATE(zvtau(jpi,jpj))
DO_2D( 0, 0, 0, 0 )
jk = mikv(ji,jj)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( pvv(ji,jj-1,jk,Kaa) + pvv(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "vtau", zvtau(:,:) )
DEALLOCATE(zvtau)
......
src/OCE/DYN/dynspg_ts.F90
View file @ 32d63bb7
......@@ -334,14 +334,14 @@ CONTAINS
! ! ------------------ !
IF( ln_bt_fw ) THEN
DO_2D( 0, 0, 0, 0 )
zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)
zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utauU(ji,jj) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtauV(ji,jj) * r1_hv(ji,jj,Kmm)
END_2D
ELSE
zztmp = r1_rho0 * r1_2
DO_2D( 0, 0, 0, 0 )
zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm)
zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utauU(ji,jj) ) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtauV(ji,jj) ) * r1_hv(ji,jj,Kmm)
END_2D
ENDIF
!
......
src/OCE/DYN/dynzdf.F90
View file @ 32d63bb7
......@@ -267,10 +267,10 @@ CONTAINS
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
#if defined key_RK3
! ! RK3: use only utau (not utau_b)
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utau(ji,jj) &
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utauU(ji,jj) &
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#else
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) &
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utauU(ji,jj) ) &
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#endif
END_2D
......@@ -397,10 +397,10 @@ CONTAINS
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
#if defined key_RK3
! ! RK3: use only vtau (not vtau_b)
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtau(ji,jj) &
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtauV(ji,jj) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#else
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtau(ji,jj) ) &
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2 * ( vtau_b(ji,jj) + vtauV(ji,jj) ) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#endif
END_2D
......
src/OCE/IOM/iom.F90
View file @ 32d63bb7
......@@ -1202,6 +1202,7 @@ CONTAINS
CHARACTER(LEN=1) :: cl_type ! local value of cd_type
LOGICAL :: ll_only3rd ! T => if kstart, kcount present then *only* use values for 3rd spatial dimension.
INTEGER :: inlev ! number of levels for 3D data
INTEGER :: ihls ! local value of the halo size
REAL(dp) :: gma, gmi
!---------------------------------------------------------------------
CHARACTER(LEN=lc) :: context
......@@ -1327,14 +1328,29 @@ CONTAINS
IF( irankpv == 1 ) ishape(1:1) = SHAPE(pv_r1d)
IF( irankpv == 2 ) ishape(1:2) = SHAPE(pv_r2d)
IF( irankpv == 3 ) ishape(1:3) = SHAPE(pv_r3d)
ix1 = 1 ; ix2 = icnt(1) ; iy1 = 1 ; iy2 = icnt(2) ! index of the array to be read
ctmp1 = 'd'
ELSE
IF( irankpv == 2 ) THEN
ishape(1:2) = SHAPE(pv_r2d(Nis0:Nie0,Njs0:Nje0 )) ; ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0)'
ENDIF
IF( irankpv == 3 ) THEN
ishape(1:3) = SHAPE(pv_r3d(Nis0:Nie0,Njs0:Nje0,:)) ; ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0,:)'
IF( irankpv == 2 ) ishape(1:2) = SHAPE(pv_r2d)
IF( irankpv == 3 ) ishape(1:3) = SHAPE(pv_r3d)
IF( ishape(1) == Ni_0 .AND. ishape(2) == Nj_0 ) THEN ! array with 0 halo
ihls = 0
ix1 = 1 ; ix2 = Ni_0 ; iy1 = 1 ; iy2 = Nj_0 ! index of the array to be read
ctmp1 = 'd(:,:'
ELSEIF( ishape(1) == jpi .AND. ishape(2) == jpj ) THEN ! array with nn_hls halos
ihls = nn_hls
ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0 ! index of the array to be read
ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0'
ELSEIF( ishape(1) == Ni_0+1 .AND. ishape(2) == Nj_0+1 ) THEN ! nn_hls = 2 and array with 1 halo
ihls = 1
ix1 = 2 ; ix2 = Ni_0+1 ; iy1 = 2 ; iy2 = Nj_0+1 ! index of the array to be read
ctmp1 = 'd(2:Ni_0+1,2:Ni_0+1'
ELSE
CALL ctl_stop( 'iom_get_123d: should have been an impossible case...' )
ENDIF
ishape(1:2) = (/ Ni_0, Nj_0 /) ! update and force ishape to match the inner domain
IF( irankpv == 3 ) ctmp1 = TRIM(ctmp1)//',:'
ctmp1 = TRIM(ctmp1)//')'
ENDIF
DO jl = 1, irankpv
WRITE( ctmp2, FMT="(', ', i1,'): ', i5,' /= icnt(', i1,'):', i5)" ) jl, ishape(jl), jl, icnt(jl)
......@@ -1347,25 +1363,22 @@ CONTAINS
!-
IF( idvar > 0 .AND. istop == nstop ) THEN ! no additional errors until this point...
!
! find the right index of the array to be read
IF( idom /= jpdom_unknown ) THEN ; ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0
ELSE ; ix1 = 1 ; ix2 = icnt(1) ; iy1 = 1 ; iy2 = icnt(2)
ENDIF
CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, pv_r3d )
IF( istop == nstop ) THEN ! no additional errors until this point...
IF(lwp) WRITE(numout,"(10x,' read ',a,' (rec: ',i6,') in ',a,' ok')") TRIM(cdvar), itime, TRIM(iom_file(kiomid)%name)
cl_type = 'T'
IF( PRESENT(cd_type) ) cl_type = cd_type
zsgn = 1._wp
IF( PRESENT(psgn ) ) zsgn = psgn
!--- overlap areas and extra hallows (mpp)
IF( PRESENT(pv_r2d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
CALL lbc_lnk( 'iom', pv_r2d, cl_type, zsgn, kfillmode = kfill )
ELSEIF( PRESENT(pv_r3d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill )
!--- halos and NP folding (NP folding to be done even if no halos)
IF( idom /= jpdom_unknown .AND. cl_type /= 'Z' .AND. ( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) ) THEN
cl_type = 'T'
IF( PRESENT(cd_type) ) cl_type = cd_type
zsgn = 1._wp
IF( PRESENT(psgn ) ) zsgn = psgn
IF( PRESENT(pv_r2d) .AND. llok ) THEN
CALL lbc_lnk( 'iom', pv_r2d, cl_type, zsgn, kfillmode = kfill, khls = ihls )
ELSEIF( PRESENT(pv_r3d) .AND. llok ) THEN
CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill, khls = ihls )
ENDIF
ENDIF
!
ELSE
......@@ -2336,14 +2349,14 @@ CONTAINS
idb(jn) = -nn_hls ! Tile data offset (halo size)
END DO
! Tile_[ij]begin are defined with respect to the processor data domain, so data_[ij]begin is added
CALL iom_set_domain_attr("grid_"//cdgrd, ntiles=nijtile, &
& tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
& tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
& tile_ni=ini(:), tile_nj=inj(:), &
& tile_data_ibegin=idb(:), tile_data_jbegin=idb(:), &
& tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
idb(:) = 0
CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", ntiles=nijtile, &
& tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
& tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
& tile_ni=ini(:), tile_nj=inj(:), &
& tile_data_ibegin=idb(:), tile_data_jbegin=idb(:), &
& tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
......@@ -2453,7 +2466,7 @@ CONTAINS
! CALL dom_ngb( -168.53_wp, 65.03_wp, ix, iy, 'T' ) ! i-line that passes through Bering Strait: Reference latitude (used in plots)
CALL dom_ngb( 180.0_wp, 90.0_wp, ix, iy, 'T' ) ! i-line that passes near the North Pole : Reference latitude (used in plots)
CALL iom_set_domain_attr("gznl", ni_glo=Ni0glo, nj_glo=Nj0glo, ibegin=mig0(Nis0)-1, jbegin=mjg0(Njs0)-1, ni=Ni_0, nj=Nj_0)
CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj)
CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin=0, data_ni=Ni_0, data_jbegin=0, data_nj=Nj_0)
CALL iom_set_domain_attr("gznl", lonvalue = real(zlon, dp), &
& latvalue = real(RESHAPE(plat(Nis0:Nie0, Njs0:Nje0),(/ Ni_0*Nj_0 /)),dp))
CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj0glo)
......
src/OCE/IOM/iom_nf90.F90
View file @ 32d63bb7
......@@ -40,6 +40,8 @@ MODULE iom_nf90
MODULE PROCEDURE iom_nf90_rp0123d_dp
END INTERFACE
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: iom_nf90.F90 14433 2021-02-11 08:06:49Z smasson $
......@@ -544,7 +546,7 @@ CONTAINS
INTEGER :: idvar ! variable id
INTEGER :: jd ! dimension loop counter
INTEGER :: ix1, ix2, iy1, iy2 ! subdomain indexes
INTEGER, DIMENSION(4) :: idimsz ! dimensions size
INTEGER, DIMENSION(3) :: ishape ! dimensions size
INTEGER, DIMENSION(4) :: idimid ! dimensions id
CHARACTER(LEN=256) :: clinfo ! info character
INTEGER :: if90id ! nf90 file identifier
......@@ -627,11 +629,9 @@ CONTAINS
itype = NF90_DOUBLE
ENDIF
IF( PRESENT(pv_r0d) ) THEN
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, &
& iom_file(kiomid)%nvid(idvar) ), clinfo )
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, iom_file(kiomid)%nvid(idvar) ), clinfo )
ELSE
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims), &
& iom_file(kiomid)%nvid(idvar) ), clinfo )
CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims), iom_file(kiomid)%nvid(idvar) ), clinfo )
ENDIF
lchunk = .false.
IF( snc4set%luse .AND. idims == 4 ) lchunk = .true.
......@@ -673,23 +673,13 @@ CONTAINS
ENDIF
! on what kind of domain must the data be written?
IF( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) THEN
idimsz(1:2) = iom_file(kiomid)%dimsz(1:2,idvar)
IF( idimsz(1) == Ni_0 .AND. idimsz(2) == Nj_0 ) THEN
ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0
ELSEIF( idimsz(1) == jpi .AND. idimsz(2) == jpj ) THEN
ix1 = 1 ; ix2 = jpi ; iy1 = 1 ; iy2 = jpj
ELSEIF( idimsz(1) == jpi .AND. idimsz(2) == jpj ) THEN
ix1 = 1 ; ix2 = jpi ; iy1 = 1 ; iy2 = jpj
ELSE
CALL ctl_stop( 'iom_nf90_rp0123d: should have been an impossible case...' )
ENDIF
! write dimension variables if it is not already done
! =============
! trick: is defined to 0 => dimension variable are defined but not yet written
IF( iom_file(kiomid)%dimsz(1, 4) == 0 ) THEN ! time_counter = 0
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 1, glamt(ix1:ix2, iy1:iy2) ), clinfo )
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 2, gphit(ix1:ix2, iy1:iy2) ), clinfo )
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 1, glamt(A2D(0)) ), clinfo )
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 2, gphit(A2D(0)) ), clinfo )
SELECT CASE (iom_file(kiomid)%comp)
CASE ('OCE')
CALL iom_nf90_check( NF90_PUT_VAR( if90id, 3, gdept_1d ), clinfo )
......@@ -704,6 +694,19 @@ CONTAINS
iom_file(kiomid)%dimsz(1, 4) = 1 ! so we don't enter this IF case any more...
IF(lwp) WRITE(numout,*) TRIM(clinfo)//' write dimension variables done'
ENDIF
IF( PRESENT(pv_r2d) ) ishape(1:2) = SHAPE(pv_r2d)
IF( PRESENT(pv_r3d) ) ishape(1:3) = SHAPE(pv_r3d)
IF( ishape(1) == Ni_0 .AND. ishape(2) == Nj_0 ) THEN ! array with 0 halo
ix1 = 1 ; ix2 = Ni_0 ; iy1 = 1 ; iy2 = Nj_0
ELSEIF( ishape(1) == jpi .AND. ishape(2) == jpj ) THEN ! array with nn_hls halos
ix1 = Nis0 ; ix2 = Nie0 ; iy1 = Njs0 ; iy2 = Nje0
ELSEIF( ishape(1) == Ni_0+1 .AND. ishape(2) == Nj_0+1 ) THEN ! nn_hls = 2 and array with 1 halo
ix1 = 2 ; ix2 = Ni_0+1 ; iy1 = 2 ; iy2 = Nj_0+1
ELSE
CALL ctl_stop( 'iom_nf90_rp0123d: should have been an impossible case...' )
ENDIF
ENDIF
! write the data
......@@ -712,7 +715,7 @@ CONTAINS
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r0d ), clinfo )
ELSEIF( PRESENT(pv_r1d) ) THEN
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r1d(:) ), clinfo )
ELSEIF( PRESENT(pv_r2d) ) THEN
ELSEIF( PRESENT(pv_r2d) ) THEN
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r2d(ix1:ix2,iy1:iy2) ), clinfo )
ELSEIF( PRESENT(pv_r3d) ) THEN
CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r3d(ix1:ix2,iy1:iy2,:) ), clinfo )
......
src/OCE/IOM/prtctl.F90
View file @ 32d63bb7
......@@ -137,9 +137,14 @@ CONTAINS
INTEGER :: jn, jl, kdir
INTEGER :: iis, iie, jjs, jje
INTEGER :: itra, inum
INTEGER, DIMENSION(4) :: ishape
REAL(2*wp) :: zsum1, zsum2, zvctl1, zvctl2
!!----------------------------------------------------------------------
!
IF( ( ktab2d_1 * ktab3d_1 * ktab4d_1 * ktab2d_2 * ktab3d_2 ) /= 0 ) THEN
CALL ctl_stop( 'prt_ctl is not working with tiles' )
ENDIF
! Arrays, scalars initialization
cl1 = ''
cl2 = ''
......@@ -157,12 +162,19 @@ CONTAINS
! Loop over each sub-domain, i.e. the total number of processors ijsplt
DO jl = 1, SIZE(nall_ictls)
! define shoter names...
iis = MAX( nall_ictls(jl), ntsi )
iie = MIN( nall_ictle(jl), ntei )
jjs = MAX( nall_jctls(jl), ntsj )
jje = MIN( nall_jctle(jl), ntej )
IF( PRESENT(tab2d_1) ) ishape(1:2) = SHAPE(tab2d_1)
IF( PRESENT(tab3d_1) ) ishape(1:3) = SHAPE(tab3d_1)
IF( PRESENT(tab4d_1) ) ishape(1:4) = SHAPE(tab4d_1)
IF( ishape(1) == jpi .AND. ishape(2) == jpj ) THEN
iis = Nis0 ; iie = Nie0 ; jjs = Njs0 ; jje = Nje0
ELSE
iis = 1 ; iie = ishape(1) ; jjs = 1 ; jje = ishape(2)
ENDIF
iis = MAX( nall_ictls(jl), iis )
iie = MIN( nall_ictle(jl), iie )
jjs = MAX( nall_jctls(jl), jjs )
jje = MIN( nall_jctle(jl), jje )
IF( PRESENT(clinfo) ) THEN ; inum = numprt_top(jl)
ELSE ; inum = numprt_oce(jl)
......@@ -188,32 +200,32 @@ CONTAINS
! 2D arrays
IF( PRESENT(tab2d_1) ) THEN
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) * mask1(iis:iie,jjs:jje,1) )
ELSE ; zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) )
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) * mask1(A2D(0),1) )
ELSE ; zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) )
ENDIF
ENDIF
IF( PRESENT(tab2d_2) ) THEN
IF( PRESENT(mask2) ) THEN ; zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) * mask2(iis:iie,jjs:jje,1) )
ELSE ; zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) )
IF( PRESENT(mask2) ) THEN ; zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) * mask2(A2D(0),1) )
ELSE ; zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) )
ENDIF
ENDIF
! 3D arrays
IF( PRESENT(tab3d_1) ) THEN
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) * mask1(iis:iie,jjs:jje,1:kdir) )
ELSE ; zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) )
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) * mask1(A2D(0),1:kdir) )
ELSE ; zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) )
ENDIF
ENDIF
IF( PRESENT(tab3d_2) ) THEN
IF( PRESENT(mask2) ) THEN ; zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) * mask2(iis:iie,jjs:jje,1:kdir) )
ELSE ; zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) )
IF( PRESENT(mask2) ) THEN ; zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) * mask2(A2D(0),1:kdir) )
ELSE ; zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) )
ENDIF
ENDIF
! 4D arrays
IF( PRESENT(tab4d_1) ) THEN
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) * mask1(iis:iie,jjs:jje,1:kdir) )
ELSE ; zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) )
IF( PRESENT(mask1) ) THEN ; zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) * mask1(A2D(0),1:kdir) )
ELSE ; zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) )
ENDIF
ENDIF
......
src/OCE/LBC/lib_mpp.F90
View file @ 32d63bb7
......@@ -143,9 +143,9 @@ MODULE lib_mpp
! Neighbourgs informations
INTEGER, PARAMETER, PUBLIC :: n_hlsmax = 3
INTEGER, DIMENSION( 8), PUBLIC :: mpinei !: 8-neighbourg MPI indexes (starting at 0, -1 if no neighbourg)
INTEGER, DIMENSION(n_hlsmax,8), PUBLIC :: mpiSnei !: 8-neighbourg Send MPI indexes (starting at 0, -1 if no neighbourg)
INTEGER, DIMENSION(n_hlsmax,8), PUBLIC :: mpiRnei !: 8-neighbourg Recv MPI indexes (starting at 0, -1 if no neighbourg)
INTEGER, DIMENSION( 8), PUBLIC :: mpinei !: 8-neighbourg MPI indexes (starting at 0, -1 if no neighbourg)
INTEGER, DIMENSION(0:n_hlsmax,8), PUBLIC :: mpiSnei !: 8-neighbourg Send MPI indexes (starting at 0, -1 if no neighbourg)
INTEGER, DIMENSION(0:n_hlsmax,8), PUBLIC :: mpiRnei !: 8-neighbourg Recv MPI indexes (starting at 0, -1 if no neighbourg)
INTEGER, PARAMETER, PUBLIC :: jpwe = 1 !: WEst
INTEGER, PARAMETER, PUBLIC :: jpea = 2 !: EAst
INTEGER, PARAMETER, PUBLIC :: jpso = 3 !: SOuth
......
src/OCE/LBC/mppini.F90
View file @ 32d63bb7
......@@ -489,7 +489,9 @@ CONTAINS
! -----------------------------------------
!
! set default neighbours
mpinei(:) = impi(:,narea)
mpinei(:) = impi(:,narea) ! should be just local but is still used in icblbc and mpp_lnk_icb_generic.h90...
mpiSnei(0,:) = -1 ! no comm if no halo (but need it for NP Folding
mpiRnei(0,:) = -1
DO jh = 1, n_hlsmax
mpiSnei(jh,:) = impi(:,narea) ! default definition
mpiRnei(jh,:) = impi(:,narea)
......@@ -1411,7 +1413,18 @@ ENDIF
Ni_0 = Nie0 - Nis0 + 1
Nj_0 = Nje0 - Njs0 + 1
!
jpkm1 = jpk-1 ! " "
jpkm1 = jpk-1
!
ntile = 0 ! Initialise "no tile" by default
nijtile = 1
ntsi = Nis0
ntsj = Njs0
ntei = Nie0
ntej = Nje0
nthl = 0
nthr = 0
nthb = 0
ntht = 0
!
END SUBROUTINE init_doloop
......