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 1069 additions and 651 deletions
src/OCE/DIA/diahsb.F90
View file @ 4dbfc9bb
......@@ -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,42 +83,73 @@ 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(17) :: zbg
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
DO ji=1,jpi
DO jj=1,jpj
ztmp(ji,jj,9 ) = - surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
ztmp(ji,jj,10) = - surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
END DO
END DO
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,9 ) = - surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
ztmp(ji,jj,10) = - surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
END_2D
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,20 +184,22 @@ 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
IF( ln_isfcav ) THEN ! ISF case
DO ji = 1, jpi
DO jj = 1, jpj
ztmp(ji,jj,12) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )
ztmp(ji,jj,13) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )
END DO
END DO
DO_2D( 0, 0, 0, 0 )
ztmp(ji,jj,12) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )
ztmp(ji,jj,13) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )
END_2D
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 +219,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,32 +357,34 @@ 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 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
END DO
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
! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
DO_3D( 0, 0, 0, 0, 1, jpk )
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_3D
frc_v = 0._wp ! volume trend due to forcing
frc_t = 0._wp ! heat content - - - -
frc_s = 0._wp ! salt content - - - -
IF( ln_linssh ) THEN
IF( ln_isfcav ) THEN
DO ji = 1, jpi
DO jj = 1, jpj
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
DO_2D( 0, 0, 0, 0 )
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_2D
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 +432,7 @@ CONTAINS
INTEGER, INTENT(in) :: Kmm ! time level index
!
INTEGER :: ierror, ios ! local integer
INTEGER :: ji, jj ! loop index
!!
NAMELIST/namhsb/ ln_diahsb
!!----------------------------------------------------------------------
......@@ -413,13 +458,13 @@ CONTAINS
! ------------------- !
! 1 - Allocate memory !
! ------------------- !
ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), &
& e3t_ini(jpi,jpj,jpk), surf(jpi,jpj), ssh_ini(jpi,jpj), tmask_ini(jpi,jpj,jpk),STAT=ierror )
ALLOCATE( hc_loc_ini(A2D(0),jpk), sc_loc_ini(A2D(0),jpk), surf_ini(A2D(0)), &
& e3t_ini(A2D(0),jpk), surf(A2D(0)), ssh_ini(A2D(0)), tmask_ini(A2D(0),jpk),STAT=ierror )
IF( ierror > 0 ) THEN
CALL ctl_stop( 'dia_hsb_init: unable to allocate hc_loc_ini' ) ; RETURN
ENDIF
IF( ln_linssh ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )
IF( ln_linssh ) ALLOCATE( ssh_hc_loc_ini(A2D(0)), ssh_sc_loc_ini(A2D(0)),STAT=ierror )
IF( ierror > 0 ) THEN
CALL ctl_stop( 'dia_hsb: unable to allocate ssh_hc_loc_ini' ) ; RETURN
ENDIF
......@@ -427,7 +472,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) * smask0_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 @ 4dbfc9bb
......@@ -54,8 +54,8 @@ CONTAINS
INTEGER :: dia_hth_alloc
!!---------------------------------------------------------------------
!
ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd26(jpi,jpj), hd28(jpi,jpj), &
& htc3(jpi,jpj), htc7(jpi,jpj), htc20(jpi,jpj), STAT=dia_hth_alloc )
ALLOCATE( hth(A2D(0)), hd20(A2D(0)), hd26(A2D(0)), hd28(A2D(0)), &
& htc3(A2D(0)), htc7(A2D(0)), htc20(A2D(0)), STAT=dia_hth_alloc )
!
CALL mpp_sum ( 'diahth', dia_hth_alloc )
IF(dia_hth_alloc /= 0) CALL ctl_stop( 'STOP', 'dia_hth_alloc: failed to allocate arrays.' )
......@@ -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)
!
......@@ -301,17 +301,20 @@ CONTAINS
!
INTEGER , INTENT(in) :: Kmm ! ocean time level index
REAL(wp), INTENT(in) :: ptem
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdept
REAL(wp), DIMENSION(A2D(0)), INTENT(out) :: pdept
!
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
!
......@@ -343,21 +346,32 @@ CONTAINS
INTEGER , INTENT(in) :: Kmm ! ocean time level index
REAL(wp), INTENT(in) :: pdep ! depth over the heat content
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pt
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc
REAL(wp), DIMENSION(A2D(0)), 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 @ 4dbfc9bb
......@@ -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 @ 4dbfc9bb
......@@ -78,7 +78,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! ocean time-step index
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in), OPTIONAL :: pvtr ! j-effective transport
REAL(wp), DIMENSION(A2D(0),jpk) , INTENT(in), OPTIONAL :: pvtr ! j-effective transport
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dia_ptr')
......@@ -110,13 +110,13 @@ CONTAINS
!!----------------------------------------------------------------------
!! ** Purpose : Calculate diagnostics and send to XIOS
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! ocean time-step index
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in), OPTIONAL :: pvtr ! j-effective transport
INTEGER , INTENT(in) :: kt ! ocean time-step index
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: pvtr ! j-effective transport (used only by PRESENT)
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
REAL(wp), DIMENSION(jpj) :: zvsum, ztsum, zssum ! 1D workspace
REAL(wp), DIMENSION(Nis0:Nie0,Njs0:Nje0) :: z2d ! 2D workspace
REAL(wp), DIMENSION( Njs0:Nje0) :: zvsum, ztsum, zssum ! 1D workspace
!
!overturning calculation
REAL(wp), DIMENSION(:,:,: ), ALLOCATABLE :: sjk, r1_sjk, v_msf ! i-mean i-k-surface and its inverse
......@@ -126,19 +126,19 @@ CONTAINS
REAL(wp), DIMENSION(:,:,: ), ALLOCATABLE :: z3dtr
!!----------------------------------------------------------------------
!
ALLOCATE( z3dtr(jpi,jpj,nbasin) )
ALLOCATE( z3dtr(Nis0:Nie0,Njs0:Nje0,nbasin) )
IF( PRESENT( pvtr ) ) THEN
IF( iom_use( 'zomsf' ) ) THEN ! effective MSF
ALLOCATE( z4d1(jpi,jpj,jpk,nbasin) )
ALLOCATE( z4d1(Nis0:Nie0,Njs0:Nje0,jpk,nbasin) )
!
DO jn = 1, nbasin ! by sub-basins
z4d1(1,:,:,jn) = pvtr_int(:,:,jp_vtr,jn) ! zonal cumulative effective transport excluding closed seas
z4d1(Nis0,:,:,jn) = pvtr_int(:,:,jp_vtr,jn) ! zonal cumulative effective transport excluding closed seas
DO jk = jpkm1, 1, -1
z4d1(1,:,jk,jn) = z4d1(1,:,jk+1,jn) - z4d1(1,:,jk,jn) ! effective j-Stream-Function (MSF)
z4d1(Nis0,:,jk,jn) = z4d1(Nis0,:,jk+1,jn) - z4d1(Nis0,:,jk,jn) ! effective j-Stream-Function (MSF)
END DO
DO ji = 2, jpi
z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
DO ji = Nis0+1, Nie0
z4d1(ji,:,:,jn) = z4d1(Nis0,:,:,jn)
ENDDO
END DO
CALL iom_put( 'zomsf', z4d1 * rc_sv )
......@@ -146,8 +146,8 @@ CONTAINS
DEALLOCATE( z4d1 )
ENDIF
IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
ALLOCATE( sjk(jpj,jpk,nbasin), r1_sjk(jpj,jpk,nbasin), v_msf(jpj,jpk,nbasin), &
& zt_jk(jpj,jpk,nbasin), zs_jk(jpj,jpk,nbasin) )
ALLOCATE( sjk( Njs0:Nje0,jpk,nbasin), r1_sjk(Njs0:Nje0,jpk,nbasin), v_msf(Njs0:Nje0,jpk,nbasin), &
& zt_jk(Njs0:Nje0,jpk,nbasin), zs_jk( Njs0:Nje0,jpk,nbasin) )
!
DO jn = 1, nbasin
sjk(:,:,jn) = pvtr_int(:,:,jp_msk,jn)
......@@ -162,16 +162,16 @@ CONTAINS
!
ENDDO
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophtove', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopstove', z3dtr )
......@@ -181,7 +181,7 @@ CONTAINS
IF( iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
! Calculate barotropic heat and salt transport here
ALLOCATE( sjk(jpj,1,nbasin), r1_sjk(jpj,1,nbasin) )
ALLOCATE( sjk(A1Dj(0),1,nbasin), r1_sjk(A1Dj(0),1,nbasin) )
!
DO jn = 1, nbasin
sjk(:,1,jn) = SUM( pvtr_int(:,:,jp_msk,jn), 2 )
......@@ -196,16 +196,16 @@ CONTAINS
!
ENDDO
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophtbtr', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopstbtr', z3dtr )
......@@ -218,28 +218,28 @@ CONTAINS
pvtr_int(:,:,:,:) = 0._wp
ELSE
IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' ) ) THEN ! i-mean i-k-surface
ALLOCATE( z4d1(jpi,jpj,jpk,nbasin), z4d2(jpi,jpj,jpk,nbasin) )
ALLOCATE( z4d1(Nis0:Nie0,Njs0:Nje0,jpk,nbasin), z4d2(Nis0:Nie0,Njs0:Nje0,jpk,nbasin) )
!
DO jn = 1, nbasin
z4d1(1,:,:,jn) = pzon_int(:,:,jp_msk,jn)
DO ji = 2, jpi
z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
z4d1(Nis0,:,:,jn) = pzon_int(:,:,jp_msk,jn)
DO ji = Nis0+1, Nie0
z4d1(ji,:,:,jn) = z4d1(Nis0,:,:,jn)
ENDDO
ENDDO
CALL iom_put( 'zosrf', z4d1 )
!
DO jn = 1, nbasin
z4d2(1,:,:,jn) = pzon_int(:,:,jp_tem,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
DO ji = 2, jpi
z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
z4d2(Nis0,:,:,jn) = pzon_int(:,:,jp_tem,jn) / MAX( z4d1(Nis0,:,:,jn), 10.e-15 )
DO ji = Nis0+1, Nie0
z4d2(ji,:,:,jn) = z4d2(Nis0,:,:,jn)
ENDDO
ENDDO
CALL iom_put( 'zotem', z4d2 )
!
DO jn = 1, nbasin
z4d2(1,:,:,jn) = pzon_int(:,:,jp_sal,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
DO ji = 2, jpi
z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
z4d2(Nis0,:,:,jn) = pzon_int(:,:,jp_sal,jn) / MAX( z4d1(Nis0,:,:,jn), 10.e-15 )
DO ji = Nis0+1, Nie0
z4d2(ji,:,:,jn) = z4d2(Nis0,:,:,jn)
ENDDO
ENDDO
CALL iom_put( 'zosal', z4d2 )
......@@ -251,16 +251,16 @@ CONTAINS
IF( iom_use( 'sophtadv' ) .OR. iom_use( 'sopstadv' ) ) THEN
!
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophtadv', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopstadv', z3dtr )
......@@ -269,16 +269,16 @@ CONTAINS
IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) THEN
!
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophtldf', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopstldf', z3dtr )
......@@ -287,16 +287,16 @@ CONTAINS
IF( iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) ) THEN
!
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophteiv', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopsteiv', z3dtr )
......@@ -304,16 +304,16 @@ CONTAINS
!
IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt ! (conversion in PW)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sophtvtr', z3dtr )
DO jn = 1, nbasin
z3dtr(1,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = 2, jpi
z3dtr(ji,:,jn) = z3dtr(1,:,jn)
z3dtr(Nis0,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram ! (conversion in Gg)
DO ji = Nis0+1, Nie0
z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
ENDDO
ENDDO
CALL iom_put( 'sopstvtr', z3dtr )
......@@ -349,8 +349,8 @@ CONTAINS
!! ** Action : pvtr_int - terms for volume streamfunction, heat/salt transport barotropic/overturning terms
!! pzon_int - terms for i mean temperature/salinity
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in), OPTIONAL :: pvtr ! j-effective transport
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(A2D(0),jpk), INTENT(in), OPTIONAL :: pvtr ! j-effective transport
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zmask ! 3D workspace
REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zts ! 4D workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: sjk, v_msf ! Zonal sum: i-k surface area, j-effective transport
......@@ -362,7 +362,7 @@ CONTAINS
IF( PRESENT( pvtr ) ) THEN
! i sum of effective j transport excluding closed seas
IF( iom_use( 'zomsf' ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
ALLOCATE( v_msf(A1Dj(nn_hls),jpk,nbasin) )
ALLOCATE( v_msf(A1Dj(0),jpk,nbasin) )
DO jn = 1, nbasin
v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) )
......@@ -374,16 +374,16 @@ CONTAINS
ENDIF
! i sum of j surface area, j surface area - temperature/salinity product on V grid
IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR. &
IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR. &
& iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
ALLOCATE( zmask(A2D(nn_hls),jpk), zts(A2D(nn_hls),jpk,jpts), &
& sjk(A1Dj(nn_hls),jpk,nbasin), &
& zt_jk(A1Dj(nn_hls),jpk,nbasin), zs_jk(A1Dj(nn_hls),jpk,nbasin) )
ALLOCATE( zmask( A2D(0),jpk ), zts( A2D(0),jpk,jpts ), &
& sjk( A1Dj(0),jpk,nbasin), &
& zt_jk(A1Dj(0),jpk,nbasin), zs_jk(A1Dj(0),jpk,nbasin) )
zmask(:,:,:) = 0._wp
zts(:,:,:,:) = 0._wp
DO_3D( 1, 1, 1, 0, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
zmask(ji,jj,jk) = vmask(ji,jj,jk) * zvfc
zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc !Tracers averaged onto V grid
......@@ -405,14 +405,14 @@ CONTAINS
ELSE
! i sum of j surface area - temperature/salinity product on T grid
IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' ) ) THEN
ALLOCATE( zmask(A2D(nn_hls),jpk), zts(A2D(nn_hls),jpk,jpts), &
& sjk(A1Dj(nn_hls),jpk,nbasin), &
& zt_jk(A1Dj(nn_hls),jpk,nbasin), zs_jk(A1Dj(nn_hls),jpk,nbasin) )
ALLOCATE( zmask( A2D(0),jpk ), zts( A2D(0),jpk,jpts ), &
& sjk( A1Dj(0),jpk,nbasin), &
& zt_jk(A1Dj(0),jpk,nbasin), zs_jk(A1Dj(0),jpk,nbasin) )
zmask(:,:,:) = 0._wp
zts(:,:,:,:) = 0._wp
DO_3D( 1, 1, 1, 1, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm)
zmask(ji,jj,jk) = tmask(ji,jj,jk) * zsfc
zts(ji,jj,jk,jp_tem) = ts(ji,jj,jk,jp_tem,Kmm) * zsfc
......@@ -434,11 +434,12 @@ CONTAINS
! i-k sum of j surface area - temperature/salinity product on V grid
IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
ALLOCATE( zts(A2D(nn_hls),jpk,jpts) )
zts(:,:,:,:) = 0._wp
DO_3D( 1, 1, 1, 0, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc !Tracers averaged onto V grid
zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
......@@ -538,11 +539,12 @@ CONTAINS
!! Wrapper for heat and salt transport calculations to calculate them for each basin
!! Called from all advection and/or diffusion routines
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ktra ! tracer index
CHARACTER(len=3) , INTENT(in) :: cptr ! transport type 'adv'/'ldf'/'eiv'
REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in) :: pvflx ! 3D input array of advection/diffusion
REAL(wp), DIMENSION(A1Dj(nn_hls),nbasin) :: zsj !
INTEGER :: jn !
INTEGER , INTENT(in) :: ktra ! tracer index
CHARACTER(len=3) , INTENT(in) :: cptr ! transport type 'adv'/'ldf'/'eiv'
! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
REAL(wp), DIMENSION(A2D(0),jpk) , INTENT(in) :: pvflx ! 3D input array of advection/diffusion
REAL(wp), DIMENSION(A1Dj(0),nbasin) :: zsj !
INTEGER :: jn !
DO jn = 1, nbasin
zsj(:,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
......@@ -576,13 +578,13 @@ CONTAINS
!!
!! ** Action : phstr
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpj,nbasin) , INTENT(inout) :: phstr !
REAL(wp), DIMENSION(A1Dj(nn_hls),nbasin), INTENT(in) :: pva !
REAL(wp), DIMENSION(Njs0:Nje0,nbasin), INTENT(inout) :: phstr !
REAL(wp), DIMENSION(A1Dj(0) ,nbasin), INTENT(in ) :: pva !
INTEGER :: jj
#if ! defined key_mpi_off
INTEGER, DIMENSION(1) :: ish1d
INTEGER, DIMENSION(2) :: ish2d
REAL(wp), DIMENSION(jpj*nbasin) :: zwork
INTEGER, DIMENSION(1) :: ish1d
INTEGER, DIMENSION(2) :: ish2d
REAL(wp), DIMENSION(:), ALLOCATABLE :: zwork
#endif
DO jj = ntsj, ntej
......@@ -591,11 +593,13 @@ CONTAINS
#if ! defined key_mpi_off
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
ish1d(1) = jpj*nbasin
ish2d(1) = jpj ; ish2d(2) = nbasin
ALLOCATE( zwork(Nj_0*nbasin) )
ish1d(1) = Nj_0*nbasin
ish2d(1) = Nj_0 ; ish2d(2) = nbasin
zwork(:) = RESHAPE( phstr(:,:), ish1d )
CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl )
phstr(:,:) = RESHAPE( zwork, ish2d )
DEALLOCATE( zwork )
ENDIF
#endif
END SUBROUTINE ptr_sum_2d
......@@ -612,13 +616,13 @@ CONTAINS
!!
!! ** Action : phstr
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpj,jpk,nbasin) , INTENT(inout) :: phstr !
REAL(wp), DIMENSION(A1Dj(nn_hls),jpk,nbasin), INTENT(in) :: pva !
INTEGER :: jj, jk
REAL(wp), DIMENSION(Njs0:Nje0,jpk,nbasin), INTENT(inout) :: phstr !
REAL(wp), DIMENSION(A1Dj(0) ,jpk,nbasin), INTENT(in ) :: pva !
INTEGER :: jj, jk
#if ! defined key_mpi_off
INTEGER, DIMENSION(1) :: ish1d
INTEGER, DIMENSION(3) :: ish3d
REAL(wp), DIMENSION(jpj*jpk*nbasin) :: zwork
INTEGER, DIMENSION(1) :: ish1d
INTEGER, DIMENSION(3) :: ish3d
REAL(wp), DIMENSION(:), ALLOCATABLE :: zwork
#endif
DO jk = 1, jpk
......@@ -629,11 +633,13 @@ CONTAINS
#if ! defined key_mpi_off
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
ish1d(1) = jpj*jpk*nbasin
ish3d(1) = jpj ; ish3d(2) = jpk ; ish3d(3) = nbasin
ALLOCATE( zwork(Nj_0*jpk*nbasin) )
ish1d(1) = Nj_0*jpk*nbasin
ish3d(1) = Nj_0 ; ish3d(2) = jpk ; ish3d(3) = nbasin
zwork(:) = RESHAPE( phstr(:,:,:), ish1d )
CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl )
phstr(:,:,:) = RESHAPE( zwork, ish3d )
DEALLOCATE( zwork )
ENDIF
#endif
END SUBROUTINE ptr_sum_3d
......@@ -651,13 +657,13 @@ CONTAINS
! nbasin has been initialized in iom_init to define the axis "basin"
!
IF( .NOT. ALLOCATED( btmsk ) ) THEN
ALLOCATE( btmsk(jpi,jpj,nbasin) , btmsk34(jpi,jpj,nbasin), &
& hstr_adv(jpj,jpts,nbasin), hstr_eiv(jpj,jpts,nbasin), &
& hstr_ove(jpj,jpts,nbasin), hstr_btr(jpj,jpts,nbasin), &
& hstr_ldf(jpj,jpts,nbasin), hstr_vtr(jpj,jpts,nbasin), STAT=ierr(1) )
ALLOCATE( btmsk(jpi,jpj,nbasin) , btmsk34(jpi,jpj,nbasin), &
& hstr_adv(Njs0:Nje0,jpts,nbasin), hstr_eiv(Njs0:Nje0,jpts,nbasin), &
& hstr_ove(Njs0:Nje0,jpts,nbasin), hstr_btr(Njs0:Nje0,jpts,nbasin), &
& hstr_ldf(Njs0:Nje0,jpts,nbasin), hstr_vtr(Njs0:Nje0,jpts,nbasin), STAT=ierr(1) )
!
ALLOCATE( pvtr_int(jpj,jpk,jpts+2,nbasin), &
& pzon_int(jpj,jpk,jpts+1,nbasin), STAT=ierr(2) )
ALLOCATE( pvtr_int(Njs0:Nje0,jpk,jpts+2,nbasin), &
& pzon_int(Njs0:Nje0,jpk,jpts+1,nbasin), STAT=ierr(2) )
!
dia_ptr_alloc = MAXVAL( ierr )
CALL mpp_sum( 'diaptr', dia_ptr_alloc )
......@@ -677,11 +683,12 @@ CONTAINS
!!
!! ** Action : - p_fval: i-k-mean poleward flux of pvflx
!!----------------------------------------------------------------------
REAL(wp), INTENT(in), DIMENSION(A2D(nn_hls),jpk) :: pvflx ! mask flux array at V-point
REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: pmsk ! Optional 2D basin mask
! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
REAL(wp), INTENT(in), DIMENSION(A2D(0),jpk) :: pvflx ! mask flux array at V-point
REAL(wp), INTENT(in), DIMENSION(jpi,jpj ) :: pmsk ! Optional 2D basin mask
!
INTEGER :: ji, jj, jk ! dummy loop arguments
REAL(wp), DIMENSION(A1Dj(nn_hls)) :: p_fval ! function value
INTEGER :: ji, jj, jk ! dummy loop arguments
REAL(wp), DIMENSION(A1Dj(0)) :: p_fval ! function value
!!--------------------------------------------------------------------
!
p_fval(:) = 0._wp
......@@ -702,11 +709,12 @@ CONTAINS
!!
!! ** Action : - p_fval: i-k-mean poleward flux of pvflx
!!----------------------------------------------------------------------
REAL(wp) , INTENT(in), DIMENSION(A2D(nn_hls)) :: pvflx ! mask flux array at V-point
REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) :: pmsk ! Optional 2D basin mask
! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
REAL(wp) , INTENT(in), DIMENSION(A2D(0) ) :: pvflx ! mask flux array at V-point
REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) :: pmsk ! Optional 2D basin mask
!
INTEGER :: ji,jj ! dummy loop arguments
REAL(wp), DIMENSION(A1Dj(nn_hls)) :: p_fval ! function value
INTEGER :: ji, jj ! dummy loop arguments
REAL(wp), DIMENSION(A1Dj(0)) :: p_fval ! function value
!!--------------------------------------------------------------------
!
p_fval(:) = 0._wp
......@@ -725,14 +733,13 @@ CONTAINS
!!
!! ** Action : - p_fval: j-cumulated sum of pva
!!----------------------------------------------------------------------
REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) :: pva ! mask flux array at V-point
REAL(wp) , INTENT(in), DIMENSION(A2D(0)) :: pva ! mask flux array at V-point
!
INTEGER :: ji,jj,jc ! dummy loop arguments
INTEGER :: ijpj ! ???
REAL(wp), DIMENSION(jpi,jpj) :: p_fval ! function value
INTEGER :: ji,jj,jc ! dummy loop arguments
INTEGER :: ijpj ! ???
REAL(wp), DIMENSION(A2D(0)) :: p_fval ! function value
!!--------------------------------------------------------------------
!
ijpj = jpj ! ???
p_fval(:,:) = 0._wp
DO jc = 1, jpnj ! looping over all processors in j axis
DO_2D( 0, 0, 0, 0 )
......@@ -756,11 +763,11 @@ CONTAINS
!!----------------------------------------------------------------------
!!
IMPLICIT none
REAL(wp) , INTENT(in), DIMENSION(A2D(nn_hls),jpk) :: pta ! mask flux array at V-point
REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) :: pmsk ! Optional 2D basin mask
REAL(wp) , INTENT(in), DIMENSION(A2D(0) ,jpk) :: pta ! mask flux array at V-point
REAL(wp) , INTENT(in), DIMENSION(jpi,jpj ) :: pmsk ! Optional 2D basin mask
!!
INTEGER :: ji, jj, jk ! dummy loop arguments
REAL(wp), DIMENSION(A1Dj(nn_hls),jpk) :: p_fval ! return function value
REAL(wp), DIMENSION(A1Dj(0),jpk) :: p_fval ! return function value
!!--------------------------------------------------------------------
!
p_fval(:,:) = 0._wp
......
src/OCE/DIA/diawri.F90
View file @ 4dbfc9bb
......@@ -122,8 +122,9 @@ CONTAINS
REAL(wp):: zztmp , zztmpx ! local scalar
REAL(wp):: zztmp2, zztmpy ! - -
REAL(wp):: ze3
REAL(wp), DIMENSION(A2D( 0)) :: z2d ! 2D workspace
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z3d ! 3D workspace
REAL(wp), DIMENSION(A2D(0)) :: z2d0 ! 2D workspace
REAL(wp), DIMENSION(A2D(0) ,jpk) :: z3d0 ! 3D workspace
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z3d ! 3D workspace
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dia_wri')
......@@ -135,8 +136,9 @@ CONTAINS
ENDIF
! initialize arrays
z2d(:,:) = 0._wp
z3d(:,:,:) = 0._wp
z2d0(:,:) = 0._wp
z3d0(:,:,:) = 0._wp
z3d (:,:,:) = 0._wp
! Output of initial vertical scale factor
CALL iom_put("e3t_0", e3t_0(:,:,:) )
......@@ -146,44 +148,44 @@ CONTAINS
!
IF ( iom_use("tpt_dep") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
z3d0(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "tpt_dep", z3d )
CALL iom_put( "tpt_dep", z3d0 )
ENDIF
! --- vertical scale factors --- !
IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN ! time-varying e3t
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3t", z3d )
IF ( iom_use("e3tdef") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = ( ( z3d(ji,jj,jk) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
END_3D
CALL iom_put( "e3tdef", z3d )
ENDIF
ENDIF
IF ( iom_use("e3u") ) THEN ! time-varying e3u
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3u(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3u" , z3d )
ENDIF
IF ( iom_use("e3v") ) THEN ! time-varying e3v
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3v(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3v" , z3d )
ENDIF
IF ( iom_use("e3w") ) THEN ! time-varying e3w
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3w(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3w" , z3d )
ENDIF
IF ( iom_use("e3f") ) THEN ! time-varying e3f caution here at Kaa
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3f(ji,jj,jk)
END_3D
CALL iom_put( "e3f" , z3d )
......@@ -213,9 +215,9 @@ CONTAINS
IF ( iom_use("sbt") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbkt(ji,jj)
z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
z2d0(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
END_2D
CALL iom_put( "sbt", z2d ) ! bottom temperature
CALL iom_put( "sbt", z2d0 ) ! bottom temperature
ENDIF
CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) ) ! 3D salinity
......@@ -223,9 +225,9 @@ CONTAINS
IF ( iom_use("sbs") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbkt(ji,jj)
z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
z2d0(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
END_2D
CALL iom_put( "sbs", z2d ) ! bottom salinity
CALL iom_put( "sbs", z2d0 ) ! bottom salinity
ENDIF
IF( .NOT.lk_SWE ) CALL iom_put( "rhop", rhop(:,:,:) ) ! 3D potential density (sigma0)
......@@ -233,16 +235,16 @@ CONTAINS
! --- momentum --- !
IF ( iom_use("taubot") ) THEN ! bottom stress
zztmp = rho0 * 0.25_wp
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_2D( 0, 0, 0, 0 )
zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * uu(ji ,jj,mbku(ji ,jj),Kmm) )**2 &
& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm) )**2 &
& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vv(ji,jj ,mbkv(ji,jj ),Kmm) )**2 &
& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm) )**2
z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
z2d0(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
!
END_2D
CALL iom_put( "taubot", z2d )
CALL iom_put( "taubot", z2d0 )
ENDIF
CALL iom_put( "uoce", uu(:,:,:,Kmm) ) ! 3D i-current
......@@ -250,9 +252,9 @@ CONTAINS
IF ( iom_use("sbu") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbku(ji,jj)
z2d(ji,jj) = uu(ji,jj,ikbot,Kmm)
z2d0(ji,jj) = uu(ji,jj,ikbot,Kmm)
END_2D
CALL iom_put( "sbu", z2d ) ! bottom i-current
CALL iom_put( "sbu", z2d0 ) ! bottom i-current
ENDIF
CALL iom_put( "voce", vv(:,:,:,Kmm) ) ! 3D j-current
......@@ -260,18 +262,15 @@ CONTAINS
IF ( iom_use("sbv") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbkv(ji,jj)
z2d(ji,jj) = vv(ji,jj,ikbot,Kmm)
z2d0(ji,jj) = vv(ji,jj,ikbot,Kmm)
END_2D
CALL iom_put( "sbv", z2d ) ! bottom j-current
CALL iom_put( "sbv", z2d0 ) ! bottom j-current
ENDIF
! ! vertical velocity
IF( ln_zad_Aimp ) THEN
IF( iom_use('woce') ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
END_3D
CALL iom_put( "woce", z3d ) ! explicit plus implicit parts
CALL iom_put( "woce", ww+wi ) ! explicit plus implicit parts
ENDIF
ELSE
CALL iom_put( "woce", ww )
......@@ -281,15 +280,15 @@ CONTAINS
! ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
IF( ln_zad_Aimp ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wi(ji,jj,jk) )
z3d0(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wi(ji,jj,jk) )
END_3D
ELSE
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ww(ji,jj,jk)
z3d0(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ww(ji,jj,jk)
END_3D
ENDIF
CALL iom_put( "w_masstr" , z3d )
IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d * z3d )
CALL iom_put( "w_masstr" , z3d0 )
IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d0 * z3d0 )
ENDIF
CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
......@@ -304,15 +303,15 @@ CONTAINS
zztmp = ts(ji,jj,1,jp_sal,Kmm)
zztmpx = (ts(ji+1,jj,1,jp_sal,Kmm) - zztmp) * r1_e1u(ji,jj) + (zztmp - ts(ji-1,jj ,1,jp_sal,Kmm)) * r1_e1u(ji-1,jj)
zztmpy = (ts(ji,jj+1,1,jp_sal,Kmm) - zztmp) * r1_e2v(ji,jj) + (zztmp - ts(ji ,jj-1,1,jp_sal,Kmm)) * r1_e2v(ji,jj-1)
z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
z2d0(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
END_2D
CALL iom_put( "sssgrad2", z2d ) ! square of module of sss gradient
CALL iom_put( "sssgrad2", z2d0 ) ! square of module of sss gradient
IF ( iom_use("sssgrad") ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = SQRT( z2d(ji,jj) )
z2d0(ji,jj) = SQRT( z2d0(ji,jj) )
END_2D
CALL iom_put( "sssgrad", z2d ) ! module of sss gradient
CALL iom_put( "sssgrad", z2d0 ) ! module of sss gradient
ENDIF
ENDIF
......@@ -321,80 +320,80 @@ CONTAINS
zztmp = ts(ji,jj,1,jp_tem,Kmm)
zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1)
z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
z2d0(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
END_2D
CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
CALL iom_put( "sstgrad2", z2d0 ) ! square of module of sst gradient
IF ( iom_use("sstgrad") ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = SQRT( z2d(ji,jj) )
z2d0(ji,jj) = SQRT( z2d0(ji,jj) )
END_2D
CALL iom_put( "sstgrad", z2d ) ! module of sst gradient
CALL iom_put( "sstgrad", z2d0 ) ! module of sst gradient
ENDIF
ENDIF
! heat and salt contents
IF( iom_use("heatc") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
END_3D
CALL iom_put( "heatc", rho0_rcp * z2d ) ! vertically integrated heat content (J/m2)
CALL iom_put( "heatc", rho0_rcp * z2d0 ) ! vertically integrated heat content (J/m2)
ENDIF
IF( iom_use("saltc") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
END_3D
CALL iom_put( "saltc", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
CALL iom_put( "saltc", rho0 * z2d0 ) ! vertically integrated salt content (PSU*kg/m2)
ENDIF
!
IF( iom_use("salt2c") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
END_3D
CALL iom_put( "salt2c", rho0 * z2d ) ! vertically integrated square of salt content (PSU2*kg/m2)
CALL iom_put( "salt2c", rho0 * z2d0 ) ! vertically integrated square of salt content (PSU2*kg/m2)
ENDIF
!
IF ( iom_use("ke") .OR. iom_use("ke_int") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
zztmpx = uu(ji-1,jj ,jk,Kmm) + uu(ji,jj,jk,Kmm)
zztmpy = vv(ji ,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm)
z3d(ji,jj,jk) = 0.25_wp * ( zztmpx*zztmpx + zztmpy*zztmpy )
z3d0(ji,jj,jk) = 0.25_wp * ( zztmpx*zztmpx + zztmpy*zztmpy )
END_3D
CALL iom_put( "ke", z3d ) ! kinetic energy
CALL iom_put( "ke", z3d0 ) ! kinetic energy
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d0(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
END_3D
CALL iom_put( "ke_int", z2d ) ! vertically integrated kinetic energy
CALL iom_put( "ke_int", z2d0 ) ! vertically integrated kinetic energy
ENDIF
!
IF ( iom_use("sKE") ) THEN ! surface kinetic energy at T point
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = 0.25_wp * ( uu(ji ,jj,1,Kmm) * uu(ji ,jj,1,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,1,Kmm) &
z2d0(ji,jj) = 0.25_wp * ( uu(ji ,jj,1,Kmm) * uu(ji ,jj,1,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,1,Kmm) &
& + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm) &
& + vv(ji,jj ,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji,jj ) * e3v(ji,jj ,1,Kmm) &
& + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * ssmask(ji,jj)
END_2D
IF ( iom_use("sKE" ) ) CALL iom_put( "sKE" , z2d )
IF ( iom_use("sKE" ) ) CALL iom_put( "sKE" , z2d0 )
ENDIF
!
IF ( iom_use("ssKEf") ) THEN ! surface kinetic energy at F point
z2d(:,:) = 0._wp ! CAUTION : only valid in SWE, not with bathymetry
z2d0(:,:) = 0._wp ! CAUTION : only valid in SWE, not with bathymetry
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = 0.25_wp * ( uu(ji,jj ,1,Kmm) * uu(ji,jj ,1,Kmm) * e1e2u(ji,jj ) * e3u(ji,jj ,1,Kmm) &
z2d0(ji,jj) = 0.25_wp * ( uu(ji,jj ,1,Kmm) * uu(ji,jj ,1,Kmm) * e1e2u(ji,jj ) * e3u(ji,jj ,1,Kmm) &
& + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm) &
& + vv(ji ,jj,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji ,jj) * e3v(ji ,jj,1,Kmm) &
& + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm) ) &
& * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj)
END_2D
CALL iom_put( "ssKEf", z2d )
CALL iom_put( "ssKEf", z2d0 )
ENDIF
!
CALL iom_put( "hdiv", hdiv ) ! Horizontal divergence
......@@ -402,31 +401,31 @@ CONTAINS
IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = rho0 * uu(ji,jj,jk,Kmm) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
z3d0(ji,jj,jk) = rho0 * uu(ji,jj,jk,Kmm) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
END_3D
CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction
CALL iom_put( "u_masstr" , z3d0 ) ! mass transport in i-direction
IF( iom_use("u_masstr_vint") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk)
z2d0(ji,jj) = z2d0(ji,jj) + z3d0(ji,jj,jk)
END_3D
CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum
CALL iom_put( "u_masstr_vint", z2d0 ) ! mass transport in i-direction vertical sum
ENDIF
IF( iom_use("u_heattr") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
zztmp = 0.5_wp * rcp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
z2d0(ji,jj) = z2d0(ji,jj) + zztmp * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
END_3D
CALL iom_put( "u_heattr", z2d ) ! heat transport in i-direction
CALL iom_put( "u_heattr", z2d0 ) ! heat transport in i-direction
ENDIF
IF( iom_use("u_salttr") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + 0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
z2d0(ji,jj) = z2d0(ji,jj) + 0.5 * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
END_3D
CALL iom_put( "u_salttr", z2d ) ! heat transport in i-direction
CALL iom_put( "u_salttr", z2d0 ) ! heat transport in i-direction
ENDIF
ENDIF
......@@ -434,41 +433,41 @@ CONTAINS
IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = rho0 * vv(ji,jj,jk,Kmm) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
z3d0(ji,jj,jk) = rho0 * vv(ji,jj,jk,Kmm) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
END_3D
CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
CALL iom_put( "v_masstr", z3d0 ) ! mass transport in j-direction
IF( iom_use("v_heattr") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
zztmp = 0.5_wp * rcp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
z2d0(ji,jj) = z2d0(ji,jj) + zztmp * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
END_3D
CALL iom_put( "v_heattr", z2d ) ! heat transport in j-direction
CALL iom_put( "v_heattr", z2d0 ) ! heat transport in j-direction
ENDIF
IF( iom_use("v_salttr") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + 0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
z2d0(ji,jj) = z2d0(ji,jj) + 0.5 * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
END_3D
CALL iom_put( "v_salttr", z2d ) ! heat transport in j-direction
CALL iom_put( "v_salttr", z2d0 ) ! heat transport in j-direction
ENDIF
ENDIF
IF( iom_use("tosmint") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
z2d0(ji,jj) = z2d0(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
END_3D
CALL iom_put( "tosmint", z2d ) ! Vertical integral of temperature
CALL iom_put( "tosmint", z2d0 ) ! Vertical integral of temperature
ENDIF
IF( iom_use("somint") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
z2d0(ji,jj) = z2d0(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
END_3D
CALL iom_put( "somint", z2d ) ! Vertical integral of salinity
CALL iom_put( "somint", z2d0 ) ! Vertical integral of salinity
ENDIF
CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
......@@ -477,45 +476,47 @@ CONTAINS
! Output of surface vorticity terms
!
CALL iom_put( "ssplavor", ff_f ) ! planetary vorticity ( f )
!
IF ( iom_use("ssrelvor") .OR. iom_use("ssEns") .OR. &
IF ( iom_use("ssplavor") .OR. iom_use("ssrelvor") .OR. iom_use("ssEns") .OR. &
& iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
!
z2d(:,:) = 0._wp
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm) &
z2d0(ji,jj) = ff_f(ji,jj)
END_2D
CALL iom_put( "ssplavor", z2d0 ) ! planetary vorticity ( f )
DO_2D( 0, 0, 0, 0 )
z2d0(ji,jj) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm) &
& - e1u(ji ,jj+1) * uu(ji ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm) ) * r1_e1e2f(ji,jj)
END_2D
CALL iom_put( "ssrelvor", z2d ) ! relative vorticity ( zeta )
CALL iom_put( "ssrelvor", z2d0 ) ! relative vorticity ( zeta )
!
IF ( iom_use("ssEns") .OR. iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
DO_2D( 0, 0, 0, 0 )
ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
ELSE ; ze3 = 0._wp
ENDIF
z2d(ji,jj) = ze3 * z2d(ji,jj)
z2d0(ji,jj) = ze3 * z2d0(ji,jj)
END_2D
CALL iom_put( "ssrelpotvor", z2d ) ! relative potential vorticity (zeta/h)
CALL iom_put( "ssrelpotvor", z2d0 ) ! relative potential vorticity (zeta/h)
!
IF ( iom_use("ssEns") .OR. iom_use("ssabspotvor") ) THEN
DO_2D( 0, 0, 0, 0 )
ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
ELSE ; ze3 = 0._wp
ENDIF
z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj)
z2d0(ji,jj) = ze3 * ff_f(ji,jj) + z2d0(ji,jj)
END_2D
CALL iom_put( "ssabspotvor", z2d ) ! absolute potential vorticity ( q )
CALL iom_put( "ssabspotvor", z2d0 ) ! absolute potential vorticity ( q )
!
IF ( iom_use("ssEns") ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = 0.5_wp * z2d(ji,jj) * z2d(ji,jj)
z2d0(ji,jj) = 0.5_wp * z2d0(ji,jj) * z2d0(ji,jj)
END_2D
CALL iom_put( "ssEns", z2d ) ! potential enstrophy ( 1/2*q2 )
CALL iom_put( "ssEns", z2d0 ) ! potential enstrophy ( 1/2*q2 )
ENDIF
ENDIF
ENDIF
......@@ -582,8 +583,8 @@ CONTAINS
INTEGER :: jn, ierror ! local integers
REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
!
REAL(wp), DIMENSION(jpi,jpj ) :: z2d ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
REAL(wp), DIMENSION(jpi,jpj ) :: z2d0 ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d0 ! 3D workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl ! ABL 3D workspace
!!----------------------------------------------------------------------
!
......@@ -868,7 +869,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 +883,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 +893,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
......@@ -937,21 +936,21 @@ CONTAINS
IF( .NOT.ln_linssh ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm)
z3d0(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm)
END_3D
CALL histwrite( nid_T, "votemper", it, z3d, ndim_T , ndex_T ) ! heat content
CALL histwrite( nid_T, "votemper", it, z3d0, ndim_T , ndex_T ) ! heat content
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm)
z3d0(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm)
END_3D
CALL histwrite( nid_T, "vosaline", it, z3d, ndim_T , ndex_T ) ! salt content
CALL histwrite( nid_T, "vosaline", it, z3d0, ndim_T , ndex_T ) ! salt content
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj ) = ts(ji,jj, 1,jp_tem,Kmm) * e3t(ji,jj, 1,Kmm)
z2d0(ji,jj ) = ts(ji,jj, 1,jp_tem,Kmm) * e3t(ji,jj, 1,Kmm)
END_2D
CALL histwrite( nid_T, "sosstsst", it, z2d, ndim_hT, ndex_hT ) ! sea surface heat content
CALL histwrite( nid_T, "sosstsst", it, z2d0, ndim_hT, ndex_hT ) ! sea surface heat content
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj ) = ts(ji,jj, 1,jp_sal,Kmm) * e3t(ji,jj, 1,Kmm)
z2d0(ji,jj ) = ts(ji,jj, 1,jp_sal,Kmm) * e3t(ji,jj, 1,Kmm)
END_2D
CALL histwrite( nid_T, "sosaline", it, z2d, ndim_hT, ndex_hT ) ! sea surface salinity content
CALL histwrite( nid_T, "sosaline", it, z2d0, ndim_hT, ndex_hT ) ! sea surface salinity content
ELSE
CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T ) ! temperature
CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T ) ! salinity
......@@ -960,41 +959,41 @@ CONTAINS
ENDIF
IF( .NOT.ln_linssh ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
z3d0(ji,jj,jk) = e3t(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
END_3D
CALL histwrite( nid_T, "vovvle3t", it, z3d , ndim_T , ndex_T ) ! level thickness
CALL histwrite( nid_T, "vovvle3t", it, z3d0 , ndim_T , ndex_T ) ! level thickness
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
z3d0(ji,jj,jk) = gdept(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
END_3D
CALL histwrite( nid_T, "vovvldep", it, z3d , ndim_T , ndex_T ) ! t-point depth
CALL histwrite( nid_T, "vovvldep", it, z3d0 , ndim_T , ndex_T ) ! t-point depth
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ( ( e3t(ji,jj,jk,Kmm) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
z3d0(ji,jj,jk) = ( ( e3t(ji,jj,jk,Kmm) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
END_3D
CALL histwrite( nid_T, "vovvldef", it, z3d , ndim_T , ndex_T ) ! level thickness deformation
CALL histwrite( nid_T, "vovvldef", it, z3d0 , ndim_T , ndex_T ) ! level thickness deformation
ENDIF
CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm) , ndim_hT, ndex_hT ) ! sea surface height
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
z2d0(ji,jj) = emp(ji,jj) - rnf(ji,jj)
END_2D
CALL histwrite( nid_T, "sowaflup", it, z2d , ndim_hT, ndex_hT ) ! upward water flux
CALL histwrite( nid_T, "sowaflup", it, z2d0 , ndim_hT, ndex_hT ) ! upward water flux
CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
! (includes virtual salt flux beneath ice
! in linear free surface case)
IF( ln_linssh ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
z2d0(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
END_2D
CALL histwrite( nid_T, "sosst_cd", it, z2d, ndim_hT, ndex_hT ) ! c/d term on sst
CALL histwrite( nid_T, "sosst_cd", it, z2d0, ndim_hT, ndex_hT ) ! c/d term on sst
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
z2d0(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
END_2D
CALL histwrite( nid_T, "sosss_cd", it, z2d, ndim_hT, ndex_hT ) ! c/d term on sss
CALL histwrite( nid_T, "sosss_cd", it, z2d0, ndim_hT, ndex_hT ) ! c/d term on sss
ENDIF
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
z2d0(ji,jj) = qsr(ji,jj) + qns(ji,jj)
END_2D
CALL histwrite( nid_T, "sohefldo", it, z2d , ndim_hT, ndex_hT ) ! total heat flux
CALL histwrite( nid_T, "sohefldo", it, z2d0 , ndim_hT, ndex_hT ) ! total heat flux
CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
IF( ALLOCATED(hmld) ) THEN ! zdf_mxl not called by SWE
CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
......@@ -1053,9 +1052,9 @@ CONTAINS
CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = erp(ji,jj) * ts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1)
z2d0(ji,jj) = erp(ji,jj) * ts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1)
END_2D
CALL histwrite( nid_T, "sosafldp", it, z2d , ndim_hT, ndex_hT ) ! salt flux damping
CALL histwrite( nid_T, "sosafldp", it, z2d0 , ndim_hT, ndex_hT ) ! salt flux damping
ENDIF
! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
......@@ -1066,18 +1065,18 @@ 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 )
z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
z3d0(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
END_3D
CALL histwrite( nid_W, "vovecrtz", it, z3d , ndim_T, ndex_T ) ! vert. current
CALL histwrite( nid_W, "vovecrtz", it, z3d0 , ndim_T, ndex_T ) ! vert. current
ELSE
CALL histwrite( nid_W, "vovecrtz", it, ww , ndim_T, ndex_T ) ! vert. current
ENDIF
......@@ -1126,8 +1125,8 @@ CONTAINS
!!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: inum
REAL(wp), DIMENSION(jpi,jpj) :: z2d
REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d
REAL(wp), DIMENSION(A2D(0)) :: z2d0
REAL(wp), DIMENSION(A2D(0),jpk) :: z3d0
!!----------------------------------------------------------------------
!
IF(lwp) THEN
......@@ -1146,9 +1145,9 @@ CONTAINS
CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm) ) ! now j-velocity
IF( ln_zad_Aimp ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
z3d0(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
END_3D
CALL iom_rstput( 0, 0, inum, 'vovecrtz', z3d ) ! now k-velocity
CALL iom_rstput( 0, 0, inum, 'vovecrtz', z3d0 ) ! now k-velocity
ELSE
CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww ) ! now k-velocity
ENDIF
......@@ -1184,26 +1183,26 @@ CONTAINS
CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
ENDIF
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
z2d0(ji,jj) = emp(ji,jj) - rnf(ji,jj)
END_2D
CALL iom_rstput( 0, 0, inum, 'sowaflup', z2d ) ! freshwater budget
CALL iom_rstput( 0, 0, inum, 'sowaflup', z2d0 ) ! freshwater budget
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
z2d0(ji,jj) = qsr(ji,jj) + qns(ji,jj)
END_2D
CALL iom_rstput( 0, 0, inum, 'sohefldo', z2d ) ! total heat flux
CALL iom_rstput( 0, 0, inum, 'sohefldo', z2d0 ) ! total heat flux
CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
IF( .NOT.ln_linssh ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
z3d0(ji,jj,jk) = gdept(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
END_3D
CALL iom_rstput( 0, 0, inum, 'vovvldep', z3d ) ! T-cell depth
CALL iom_rstput( 0, 0, inum, 'vovvldep', z3d0 ) ! T-cell depth
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
z3d0(ji,jj,jk) = e3t(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
END_3D
CALL iom_rstput( 0, 0, inum, 'vovvle3t', z3d ) ! T-cell thickness
CALL iom_rstput( 0, 0, inum, 'vovvle3t', z3d0 ) ! T-cell thickness
END IF
IF( ln_wave .AND. ln_sdw ) THEN
CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
......
src/OCE/DIU/diu_bulk.F90
View file @ 4dbfc9bb
......@@ -64,7 +64,7 @@ CONTAINS
IF( ln_diurnal ) THEN
!
ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) )
ALLOCATE( x_dsst(A2D(0)), x_solfrac(A2D(0)) )
!
x_solfrac = 0._wp ! Initialise the solar fraction
x_dsst = 0._wp
......@@ -92,25 +92,25 @@ CONTAINS
!! ** Reference : Refinements to a prognostic scheme of skin sea surface
!! temperature, Takaya et al, JGR, 2010
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! time step
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: psolflux ! solar flux (Watts)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: pqflux ! heat (non-solar) flux (Watts)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: ptauflux ! wind stress (kg/ m s^2)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: prho ! water density (kg/m^3)
REAL(wp) , INTENT(in) :: p_rdt ! time-step
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pLa ! Langmuir number
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pthick ! warm layer thickness (m)
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pcoolthick ! cool skin thickness (m)
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pmu ! mu parameter
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: p_hflux_bkginc ! increment to the heat flux
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: p_fvel_bkginc ! increment to the friction velocity
INTEGER , INTENT(in) :: kt ! time step
REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: psolflux ! solar flux (Watts)
REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: pqflux ! heat (non-solar) flux (Watts)
REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: ptauflux ! wind stress (kg/ m s^2)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: prho ! water density (kg/m^3)
REAL(wp) , INTENT(in) :: p_rdt ! time-step
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pLa ! Langmuir number
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pthick ! warm layer thickness (m)
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pcoolthick ! cool skin thickness (m)
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pmu ! mu parameter
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: p_hflux_bkginc ! increment to the heat flux
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: p_fvel_bkginc ! increment to the friction velocity
!
INTEGER :: ji,jj
LOGICAL :: ll_calcfrac
REAL(wp), DIMENSION(jpi,jpj) :: z_fvel ! friction velocity
REAL(wp), DIMENSION(jpi,jpj) :: zthick, zcoolthick, zmu, zla
REAL(wp), DIMENSION(jpi,jpj) :: z_abflux ! absorbed flux
REAL(wp), DIMENSION(jpi,jpj) :: z_fla ! Langmuir function value
REAL(wp), DIMENSION(A2D(0)) :: z_fvel ! friction velocity
REAL(wp), DIMENSION(A2D(0)) :: zthick, zcoolthick, zmu, zla
REAL(wp), DIMENSION(A2D(0)) :: z_abflux ! absorbed flux
REAL(wp), DIMENSION(A2D(0)) :: z_fla ! Langmuir function value
!!----------------------------------------------------------------------
! Set optional arguments to their defaults
......@@ -129,14 +129,14 @@ CONTAINS
! If not done already, calculate the solar fraction
IF ( kt==nit000 ) THEN
DO_2D( 1, 1, 1, 1 )
IF( ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) &
DO_2D( 0, 0, 0, 0 )
IF( ( x_solfrac(ji,jj) == 0._wp ) .AND. ( smask0(ji,jj) == 1._wp ) ) &
& x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
END_2D
ENDIF
! convert solar flux and heat flux to absorbed flux
WHERE ( tmask(:,:,1) == 1._wp)
WHERE ( smask0(:,:) == 1._wp)
z_abflux(:,:) = ( x_solfrac(:,:) * psolflux (:,:)) + pqflux(:,:)
ELSEWHERE
z_abflux(:,:) = 0._wp
......@@ -147,8 +147,8 @@ CONTAINS
ENDWHERE
! Calculate the friction velocity
WHERE ( (ptauflux /= 0) .AND. ( tmask(:,:,1) == 1.) )
z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(:,:) )
WHERE ( (ptauflux(:,:) /= 0) .AND. ( smask0(:,:) == 1.) )
z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(A2D(0)) )
ELSEWHERE
z_fvel(:,:) = 0._wp
ENDWHERE
......@@ -157,7 +157,7 @@ CONTAINS
! Calculate the Langmuir function value
WHERE ( tmask(:,:,1) == 1.)
WHERE ( smask0(:,:) == 1.)
z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) )
ELSEWHERE
z_fla(:,:) = 0._wp
......@@ -176,16 +176,16 @@ CONTAINS
IMPLICIT NONE
! Function definition
REAL(wp), DIMENSION(jpi,jpj) :: t_imp
REAL(wp), DIMENSION(A2D(0)) :: t_imp
! Dummy variables
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst ! Delta SST
REAL(wp), INTENT(IN) :: p_rdt ! Time-step
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux ! Heat forcing
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel ! Friction velocity
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fla ! Langmuir number
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pmu ! Structure parameter
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pthick ! Layer thickness
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: prho ! Water density
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_dsst ! Delta SST
REAL(wp), INTENT(IN) :: p_rdt ! Time-step
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_abflux ! Heat forcing
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_fvel ! Friction velocity
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_fla ! Langmuir number
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: pmu ! Structure parameter
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: pthick ! Layer thickness
REAL(wp), DIMENSION(jpi,jpj),INTENT(IN) :: prho ! Water density
! Local variables
REAL(wp) :: z_olength ! Obukhov length
......@@ -198,10 +198,10 @@ CONTAINS
INTEGER :: ji,jj
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
! Only calculate outside tmask
IF ( tmask(ji,jj,1) /= 1._wp ) THEN
IF ( smask0(ji,jj) /= 1._wp ) THEN
t_imp(ji,jj) = 0._wp
CYCLE
END IF
......
src/OCE/DIU/diu_coolskin.F90
View file @ 4dbfc9bb
......@@ -59,7 +59,7 @@ MODULE diu_coolskin
!! ** Reference :
!!
!!----------------------------------------------------------------------
ALLOCATE( x_csdsst(jpi,jpj), x_csthick(jpi,jpj) )
ALLOCATE( x_csdsst(A2D(0)), x_csthick(A2D(0)) )
x_csdsst = 0.
x_csthick = 0.
!
......@@ -77,16 +77,16 @@ MODULE diu_coolskin
!! ** Reference :
!!----------------------------------------------------------------------
! Dummy variables
REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psqflux ! Heat (non-solar)(Watts)
REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pstauflux ! Wind stress (kg/ m s^2)
REAL(wp), INTENT(IN), DIMENSION(A2D(0)) :: psqflux ! Heat (non-solar)(Watts)
REAL(wp), INTENT(IN), DIMENSION(A2D(0)) :: pstauflux ! Wind stress (kg/ m s^2)
REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psrho ! Water density (kg/m^3)
REAL(wp), INTENT(IN) :: pDt ! Time-step
! Local variables
REAL(wp), DIMENSION(jpi,jpj) :: z_fv ! Friction velocity
REAL(wp), DIMENSION(jpi,jpj) :: z_gamma ! Dimensionless function of wind speed
REAL(wp), DIMENSION(jpi,jpj) :: z_lamda ! Sauders (dimensionless) proportionality constant
REAL(wp), DIMENSION(jpi,jpj) :: z_wspd ! Wind speed (m/s)
REAL(wp), DIMENSION(A2D(0)) :: z_fv ! Friction velocity
REAL(wp), DIMENSION(A2D(0)) :: z_gamma ! Dimensionless function of wind speed
REAL(wp), DIMENSION(A2D(0)) :: z_lamda ! Sauders (dimensionless) proportionality constant
REAL(wp), DIMENSION(A2D(0)) :: z_wspd ! Wind speed (m/s)
REAL(wp) :: z_ztx ! Temporary u wind stress
REAL(wp) :: z_zty ! Temporary v wind stress
REAL(wp) :: z_zmod ! Temporary total wind stress
......@@ -96,10 +96,10 @@ MODULE diu_coolskin
!
IF( .NOT. (ln_blk .OR. ln_abl) ) CALL ctl_stop("diu_coolskin.f90: diurnal flux processing only implemented for bulk forcing")
!
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
!
! Calcualte wind speed from wind stress and friction velocity
IF( tmask(ji,jj,1) == 1. .AND. pstauflux(ji,jj) /= 0 .AND. psrho(ji,jj) /=0 ) THEN
IF( smask0(ji,jj) == 1. .AND. pstauflux(ji,jj) /= 0 .AND. psrho(ji,jj) /=0 ) THEN
z_fv(ji,jj) = SQRT( pstauflux(ji,jj) / psrho(ji,jj) )
z_wspd(ji,jj) = SQRT( pstauflux(ji,jj) / ( pp_cda * pp_rhoa ) )
ELSE
......@@ -108,28 +108,28 @@ MODULE diu_coolskin
ENDIF
!
! Calculate gamma function which is dependent upon wind speed
IF( tmask(ji,jj,1) == 1. ) THEN
IF( smask0(ji,jj) == 1. ) THEN
IF( ( z_wspd(ji,jj) <= 7.5 ) ) z_gamma(ji,jj) = ( 0.2 * z_wspd(ji,jj) ) + 0.5
IF( ( z_wspd(ji,jj) > 7.5 ) .AND. ( z_wspd(ji,jj) < 10. ) ) z_gamma(ji,jj) = ( 1.6 * z_wspd(ji,jj) ) - 10.
IF( ( z_wspd(ji,jj) >= 10. ) ) z_gamma(ji,jj) = 6.
ENDIF
!
! Calculate lamda function
IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 ) THEN
IF( smask0(ji,jj) == 1. .AND. z_fv(ji,jj) /= 0 ) THEN
z_lamda(ji,jj) = ( z_fv(ji,jj) * pp_k * pp_C ) / ( z_gamma(ji,jj) * psrho(ji,jj) * pp_cw * pp_h * pp_v )
ELSE
z_lamda(ji,jj) = 0.
ENDIF
!
! Calculate the cool skin thickness - only when heat flux is out of the ocean
IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 .AND. psqflux(ji,jj) < 0 ) THEN
IF( smask0(ji,jj) == 1. .AND. z_fv(ji,jj) /= 0 .AND. psqflux(ji,jj) < 0 ) THEN
x_csthick(ji,jj) = ( z_lamda(ji,jj) * pp_v ) / z_fv(ji,jj)
ELSE
x_csthick(ji,jj) = 0.
ENDIF
!
! Calculate the cool skin correction - only when the heat flux is out of the ocean
IF( tmask(ji,jj,1) == 1. .AND. x_csthick(ji,jj) /= 0. .AND. psqflux(ji,jj) < 0. ) THEN
IF( smask0(ji,jj) == 1. .AND. x_csthick(ji,jj) /= 0. .AND. psqflux(ji,jj) < 0. ) THEN
x_csdsst(ji,jj) = ( psqflux(ji,jj) * x_csthick(ji,jj) ) / pp_k
ELSE
x_csdsst(ji,jj) = 0.
......
src/OCE/DIU/step_diu.F90
View file @ 4dbfc9bb
......@@ -46,8 +46,7 @@ MODULE step_diu
!!----------------------------------------------------------------------
INTEGER :: jk ! dummy loop indices
INTEGER :: indic ! error indicator if < 0
REAL(wp), DIMENSION(jpi,jpj) :: z_fvel_bkginc, z_hflux_bkginc
INTEGER :: Nbb, Nnn, Naa, Nrhs ! local definitions as placeholders for now
INTEGER :: Nbb, Nnn, Naa, Nrhs ! local definitions as placeholders for now
!! ---------------------------------------------------------------------
IF(ln_diurnal_only) THEN
......
src/OCE/DOM/dom_oce.F90
View file @ 4dbfc9bb
......@@ -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
!! ----------------------------
......@@ -74,10 +76,10 @@ MODULE dom_oce
INTEGER :: nn_ltile_i, nn_ltile_j
! Domain tiling
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntsi_a !: start of internal part of tile domain
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntsj_a !
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntei_a !: end of internal part of tile domain
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntej_a !
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntsi_a !: start of internal part of tile domain
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntsj_a !
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntei_a !: end of internal part of tile domain
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntej_a !
LOGICAL, PUBLIC :: l_istiled ! whether tiling is currently active or not
! !: domain MPP decomposition parameters
......@@ -85,32 +87,30 @@ MODULE dom_oce
INTEGER , PUBLIC :: narea !: number for local area (starting at 1) = MPI rank + 1
INTEGER, PUBLIC :: nidom !: IOIPSL things...
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mig !: local ==> global domain, including halos (jpiglo), i-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mjg !: local ==> global domain, including halos (jpjglo), j-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mig0 !: local ==> global domain, excluding halos (Ni0glo), i-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mjg0 !: local ==> global domain, excluding halos (Nj0glo), j-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mi0, mi1 !: global, including halos (jpiglo) ==> local domain i-index
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mig !: local ==> global domain, i-index
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mjg !: local ==> global domain, j-index
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mi0, mi1 !: global ==> local domain, i-index
! !: (mi0=1 and mi1=0 if global index not in local domain)
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mj0, mj1 !: global, including halos (jpjglo) ==> local domain j-index
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mj0, mj1 !: global ==> local domain, j-index
! !: (mj0=1 and mj1=0 if global index not in local domain)
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: nfimpp, nfproc, nfjpi
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: nfimpp, nfproc, nfjpi, nfni_0
!!----------------------------------------------------------------------
!! horizontal curvilinear coordinate and scale factors
!! ---------------------------------------------------------------------
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: glamt , glamu, glamv , glamf !: longitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: gphit , gphiu, gphiv , gphif !: latitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t , e2t , r1_e1t, r1_e2t !: t-point horizontal scale factors [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u , e2u , r1_e1u, r1_e2u !: horizontal scale factors at u-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v , e2v , r1_e1v, r1_e2v !: horizontal scale factors at v-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f , e2f , r1_e1f, r1_e2f !: horizontal scale factors at f-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: glamt , glamu, glamv , glamf !: longitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: gphit , gphiu, gphiv , gphif !: latitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1t , e2t , r1_e1t, r1_e2t !: t-point horizontal scale factors [m]
REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1u , e2u , r1_e1u, r1_e2u !: horizontal scale factors at u-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1v , e2v , r1_e1v, r1_e2v !: horizontal scale factors at v-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1f , e2f , r1_e1f, r1_e2f !: horizontal scale factors at f-point [m]
!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2t , r1_e1e2t !: associated metrics at t-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2u , r1_e1e2u , e2_e1u !: associated metrics at u-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2v , r1_e1e2v , e1_e2v !: associated metrics at v-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2f , r1_e1e2f !: associated metrics at f-point
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2t , r1_e1e2t !: associated metrics at t-point
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2u , r1_e1e2u , e2_e1u !: associated metrics at u-point
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2v , r1_e1e2v , e1_e2v !: associated metrics at v-point
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2f , r1_e1e2f !: associated metrics at f-point
!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff_f , ff_t !: Coriolis factor at f- & t-points [1/s]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ff_f , ff_t !: Coriolis factor at f- & t-points [1/s]
!!----------------------------------------------------------------------
!! vertical coordinate and scale factors
......@@ -130,75 +130,77 @@ MODULE dom_oce
LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate
LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF
! ! reference scale factors
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0 !: t- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_0 !: u- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v_0 !: v- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f_0 !: f- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_0 !: w- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_0 !: uw-vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_0 !: vw-vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3t_0 !: t- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3u_0 !: u- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3v_0 !: v- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3f_0 !: f- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3w_0 !: w- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3uw_0 !: uw-vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3vw_0 !: vw-vert. scale factor [m]
! ! time-dependent scale factors (domvvl)
#if defined key_qco || defined key_linssh
#else
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e3t, e3u, e3v, e3w, e3uw, e3vw !: vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f !: F-point vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:) :: e3t, e3u, e3v, e3w, e3uw, e3vw !: vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3f !: F-point vert. scale factor [m]
#endif
! ! time-dependent ratio ssh / h_0 (domqco)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: r3t, r3u, r3v !: time-dependent ratio at t-, u- and v-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3f !: mid-time-level ratio at f-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3t_f, r3u_f, r3v_f !: now time-filtered ratio at t-, u- and v-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: r3t, r3u, r3v !: time-dependent ratio at t-, u- and v-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: r3f !: mid-time-level ratio at f-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: r3t_f, r3u_f, r3v_f !: now time-filtered ratio at t-, u- and v-point [-]
! ! reference depths of cells
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0 !: t- depth [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdepw_0 !: w- depth [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: gdept_0 !: t- depth [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: gdepw_0 !: w- depth [m]
! ! time-dependent depths of cells (domvvl)
#if defined key_qco || defined key_linssh
#else
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w_0, gde3w !: w- depth (sum of e3w) [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: gdept, gdepw
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: gde3w_0, gde3w !: w- depth (sum of e3w) [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:) :: gdept, gdepw
#endif
! ! reference heights of ocean water column and its inverse
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0, r1_ht_0 !: t-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0, r1_hu_0 !: u-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0, r1_hv_0 !: v-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hf_0, r1_hf_0 !: f-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht_0, r1_ht_0 !: t-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hu_0, r1_hu_0 !: u-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hv_0, r1_hv_0 !: v-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hf_0, r1_hf_0 !: f-depth [m] and [1/m]
! ! time-dependent heights of ocean water column (domvvl)
#if defined key_qco || defined key_linssh
#else
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht !: t-points [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hu, r1_hu !: u-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hv, r1_hv !: v-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht !: t-points [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: hu, r1_hu !: u-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: hv, r1_hv !: v-depth [m] and [1/m]
#endif
INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1)
INTEGER, PUBLIC :: nlb10 !: shallowest W level Bellow ~10m (nla10 + 1)
!! 1D reference vertical coordinate
!! =-----------------====------
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: e3t_1d , e3w_1d !: reference vertical scale factors at T- and W-pts (m)
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m)
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: e3t_1d , e3w_1d !: reference vertical scale factors at T- and W-pts (m)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfdep, bathy
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: risfdep, bathy
!!----------------------------------------------------------------------
!! masks, top and bottom ocean point position
!! ---------------------------------------------------------------------
!!gm Proposition of new name for top/bottom vertical indices
! INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mtk_t, mtk_u, mtk_v !: top first wet T-, U-, and V-level (ISF)
! INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbk_t, mbk_u, mbk_v !: bottom last wet T-, U-, and V-level
! INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mtk_t, mtk_u, mtk_v !: top first wet T-, U-, and V-level (ISF)
! INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mbk_t, mbk_u, mbk_v !: bottom last wet T-, U-, and V-level
!!gm
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt, mbku, mbkv, mbkf !: bottom last wet T-, U-, V- and F-level
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior (excluding halos+duplicated points) domain T-point mask
INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mbkt, mbku, mbkv, mbkf !: bottom last wet T-, U-, V- and F-level
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tmask_i !: interior (excluding halos+duplicated points) domain T-point mask
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: top first wet T-, U-, V-, F-level (ISF)
INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: top first wet T-, U-, V-, F-level (ISF)
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
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: fe3mask !: land/ocean mask at F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_upd, umask_upd, vmask_upd !: land/ocean mask at F-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: smask0 !: surface mask at T-pts on inner domain
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: smask0_i !: equivalent of tmask_i for inner domain
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, wmask, fmask !: land/ocean mask at T-, U-, V-, W- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:), TARGET :: wumask, wvmask !: land/ocean mask at WU- and WV-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:), TARGET :: fe3mask !: land/ocean mask at F-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tmask_upd, umask_upd, vmask_upd !: land/ocean mask at F-pts
!!----------------------------------------------------------------------
!! calendar variables
......@@ -326,7 +328,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(Nis0-(0):Nie0+(0),Njs0-(0):Nje0+(0)) , smask0_i(Nis0-(0):Nie0+(0),Njs0-(0):Nje0+(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 @ 4dbfc9bb
......@@ -144,7 +144,8 @@ CONTAINS
tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
END_3D
ENDIF
smask0(:,:) = 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
......@@ -194,7 +195,8 @@ CONTAINS
! --------------------
!
CALL dom_uniq( tmask_i, 'T' )
tmask_i(:,:) = ssmask(:,:) * tmask_i(:,:)
tmask_i (:,:) = ssmask(:,:) * tmask_i(:,:)
smask0_i(:,:) = tmask_i(A2D(0))
! Lateral boundary conditions on velocity (modify fmask)
! ---------------------------------------
......
src/OCE/DOM/domtile.F90
View file @ 4dbfc9bb
......@@ -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/DOM/domvvl.F90
View file @ 4dbfc9bb
......@@ -282,8 +282,8 @@ CONTAINS
IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
ii0 = 103 + nn_hls - 1 ; ii1 = 111 + nn_hls - 1
ij0 = 128 + nn_hls ; ij1 = 135 + nn_hls
frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp
frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rn_Dt
frq_rst_e3t( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) ) = 0.0_wp
frq_rst_hdv( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) ) = 1.e0_wp / rn_Dt
ENDIF
ENDIF
ENDIF
......
src/OCE/DOM/domzgr.F90
View file @ 4dbfc9bb
......@@ -130,14 +130,14 @@ CONTAINS
!
zmsk(:,:) = 1._wp ! default: no closed boundaries
IF( .NOT. l_Iperio ) THEN ! E-W closed:
zmsk( mi0( 1+nn_hls):mi1( 1+nn_hls),:) = 0._wp ! first column of inner global domain at 0
zmsk( mi0(jpiglo-nn_hls):mi1(jpiglo-nn_hls),:) = 0._wp ! last column of inner global domain at 0
zmsk( mi0( 1+nn_hls,nn_hls):mi1( 1+nn_hls,nn_hls),:) = 0._wp ! first column of inner global domain at 0
zmsk( mi0(jpiglo-nn_hls,nn_hls):mi1(jpiglo-nn_hls,nn_hls),:) = 0._wp ! last column of inner global domain at 0
ENDIF
IF( .NOT. l_Jperio ) THEN ! S closed:
zmsk(:,mj0( 1+nn_hls):mj1( 1+nn_hls) ) = 0._wp ! first line of inner global domain at 0
zmsk(:,mj0( 1+nn_hls,nn_hls):mj1( 1+nn_hls,nn_hls) ) = 0._wp ! first line of inner global domain at 0
ENDIF
IF( .NOT. ( l_Jperio .OR. l_NFold ) ) THEN ! N closed:
zmsk(:,mj0(jpjglo-nn_hls):mj1(jpjglo-nn_hls) ) = 0._wp ! last line of inner global domain at 0
zmsk(:,mj0(jpjglo-nn_hls,nn_hls):mj1(jpjglo-nn_hls,nn_hls) ) = 0._wp ! last line of inner global domain at 0
ENDIF
CALL lbc_lnk( 'usrdef_zgr', zmsk, 'T', 1. ) ! set halos
k_top(:,:) = k_top(:,:) * NINT( zmsk(:,:) )
......
src/OCE/DOM/dtatsd.F90
View file @ 4dbfc9bb
......@@ -161,8 +161,8 @@ CONTAINS
ij0 = 101 + nn_hls ; ij1 = 109 + nn_hls ! Reduced T & S in the Alboran Sea
ii0 = 141 + nn_hls - 1 ; ii1 = 155 + nn_hls - 1
IF( sf_tsd(jp_tem)%ln_tint .OR. irec_n(jp_tem) /= irec_b(jp_tem) ) THEN
DO jj = mj0(ij0), mj1(ij1)
DO ji = mi0(ii0), mi1(ii1)
DO jj = mj0(ij0,nn_hls), mj1(ij1,nn_hls)
DO ji = mi0(ii0,nn_hls), mi1(ii1,nn_hls)
sf_tsd(jp_tem)%fnow(ji,jj,13:13) = sf_tsd(jp_tem)%fnow(ji,jj,13:13) - 0.20_wp
sf_tsd(jp_tem)%fnow(ji,jj,14:15) = sf_tsd(jp_tem)%fnow(ji,jj,14:15) - 0.35_wp
sf_tsd(jp_tem)%fnow(ji,jj,16:25) = sf_tsd(jp_tem)%fnow(ji,jj,16:25) - 0.40_wp
......@@ -172,8 +172,8 @@ CONTAINS
ENDIF
!
IF( sf_tsd(jp_sal)%ln_tint .OR. irec_n(jp_sal) /= irec_b(jp_sal) ) THEN
DO jj = mj0(ij0), mj1(ij1)
DO ji = mi0(ii0), mi1(ii1)
DO jj = mj0(ij0,nn_hls), mj1(ij1,nn_hls)
DO ji = mi0(ii0,nn_hls), mi1(ii1,nn_hls)
sf_tsd(jp_sal)%fnow(ji,jj,13:13) = sf_tsd(jp_sal)%fnow(ji,jj,13:13) - 0.15_wp
sf_tsd(jp_sal)%fnow(ji,jj,14:15) = sf_tsd(jp_sal)%fnow(ji,jj,14:15) - 0.25_wp
sf_tsd(jp_sal)%fnow(ji,jj,16:17) = sf_tsd(jp_sal)%fnow(ji,jj,16:17) - 0.30_wp
......@@ -185,9 +185,9 @@ CONTAINS
!
ij0 = 87 + nn_hls ; ij1 = 96 + nn_hls ! Reduced temperature in Red Sea
ii0 = 148 + nn_hls - 1 ; ii1 = 160 + nn_hls - 1
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 4:10 ) = 7.0_wp
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp
sf_tsd(jp_tem)%fnow( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 4:10 ) = 7.0_wp
sf_tsd(jp_tem)%fnow( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 11:13 ) = 6.5_wp
sf_tsd(jp_tem)%fnow( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 14:20 ) = 6.0_wp
ENDIF
ENDIF
!!gm end
......
src/OCE/DYN/dynadv.F90
View file @ 4dbfc9bb
......@@ -7,6 +7,7 @@ MODULE dynadv
!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
!! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
!! 4.0 ! 2017-07 (G. Madec) add a linear dynamics option
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -50,7 +51,7 @@ MODULE dynadv
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!---------------------------------------------------------------------
!! *** ROUTINE dyn_adv ***
!!
......@@ -64,7 +65,6 @@ CONTAINS
!! (see dynvor module).
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time step and level indices
INTEGER , OPTIONAL , INTENT(in ) :: no_zad ! no vertical advection compotation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum Eq.
!!----------------------------------------------------------------------
......@@ -73,12 +73,23 @@ CONTAINS
!
SELECT CASE( n_dynadv ) !== compute advection trend and add it to general trend ==!
CASE( np_VEC_c2 ) != vector form =!
CALL dyn_keg ( kt, nn_dynkeg , Kmm, puu, pvv, Krhs ) ! horizontal gradient of kinetic energy
CALL dyn_zad ( kt , Kmm, puu, pvv, Krhs ) ! vertical advection
! !* horizontal gradient of kinetic energy
IF (nn_hls==1) THEN ! halo 1 case
CALL dyn_keg_hls1( kt, nn_dynkeg , Kmm, puu, pvv, Krhs ) ! lbc needed with Hollingsworth scheme
ELSE ! halo 2 case
CALL dyn_keg ( kt, nn_dynkeg , Kmm, puu, pvv, Krhs )
ENDIF
CALL dyn_zad ( kt , Kmm, puu, pvv, Krhs ) !* vertical advection
!
CASE( np_FLX_c2 ) != flux form =!
CALL dyn_adv_cen2( kt , Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) ! 2nd order centered scheme
CASE( np_FLX_ubs )
CALL dyn_adv_ubs ( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) ! 3rd order UBS scheme (UP3)
CALL dyn_adv_cen2( kt , Kmm, puu, pvv, Krhs, pau, pav, paw ) !* 2nd order centered scheme
!
CASE( np_FLX_ubs ) !* 3rd order UBS scheme (UP3)
IF (nn_hls==1) THEN ! halo 1 case
CALL dyn_adv_ubs_hls1( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
ELSE ! halo 2 case
CALL dyn_adv_ubs ( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
ENDIF
END SELECT
!
IF( ln_timing ) CALL timing_stop( 'dyn_adv' )
......
src/OCE/DYN/dynadv_cen2.F90
View file @ 4dbfc9bb
......@@ -6,6 +6,7 @@ MODULE dynadv_cen2
!!======================================================================
!! History : 2.0 ! 2006-08 (G. Madec, S. Theetten) Original code
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -35,7 +36,7 @@ MODULE dynadv_cen2
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_cen2 ***
!!
......@@ -51,15 +52,17 @@ CONTAINS
!! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kmm, Krhs ! ocean time-step and level indices
INTEGER , OPTIONAL , INTENT(in ) :: no_zad ! no vertical advection computation
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zzu, zzv ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t, zfu_f, zfu_uw, zfu
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t, zfv_f, zfv_vw, zfv, zfw
REAL(wp) :: zzu, zzfu_kp1 ! local scalars
REAL(wp) :: zzv, zzfv_kp1 ! - -
REAL(wp), DIMENSION(A2D(1)) :: zfu_t, zfu_f, zfu
REAL(wp), DIMENSION(A2D(1)) :: zfv_t, zfv_f, zfv
REAL(wp), DIMENSION(A2D(1)) :: zfu_uw, zfv_vw, zfw
REAL(wp), DIMENSION(:,:,:) , POINTER :: zpt_u, zpt_v, zpt_w
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zu_trd, zv_trd
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......@@ -71,8 +74,9 @@ CONTAINS
ENDIF
!
IF( l_trddyn ) THEN ! trends: store the input trends
zfu_uw(:,:,:) = puu(:,:,:,Krhs)
zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
ALLOCATE( zu_trd(A2D(0),jpkm1), zv_trd(A2D(0),jpkm1) )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
......@@ -89,84 +93,86 @@ CONTAINS
!
DO jk = 1, jpkm1 ! horizontal transport
DO_2D( 1, 1, 1, 1 )
zfu(ji,jj,jk) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
zfu(ji,jj) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
zfv(ji,jj) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
END_2D
DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes (at T- and F-point)
zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
zfv_f(ji ,jj ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) )
zfu_f(ji ,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) )
zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
zfu_t(ji+1,jj ) = ( zfu(ji,jj) + zfu(ji+1,jj) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
zfv_f(ji ,jj ) = ( zfv(ji,jj) + zfv(ji+1,jj) ) * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) )
zfu_f(ji ,jj ) = ( zfu(ji,jj) + zfu(ji,jj+1) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) )
zfv_t(ji ,jj+1) = ( zfv(ji,jj) + zfv(ji,jj+1) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
END_2D
DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) &
& + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) &
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj) - zfu_t(ji,jj ) &
& + zfv_f(ji ,jj) - zfv_f(ji,jj-1) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) &
& + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) &
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ) - zfu_f(ji-1,jj) &
& + zfv_t(ji,jj+1) - zfv_t(ji ,jj) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
END DO
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:)
zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:)
CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt, Kmm )
zfu_t(:,:,:) = puu(:,:,:,Krhs)
zfv_t(:,:,:) = pvv(:,:,:,Krhs)
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
IF( PRESENT( no_zad ) ) THEN !== No vertical advection ==! (except if linear free surface)
! ==
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,1,Kmm)
pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,1,Kmm)
END_2D
ENDIF
!
ELSE !== Vertical advection ==!
! !== Vertical advection ==!
!
! ! surface vertical fluxes
!
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
ELSE ! non linear free: surface advective fluxes set to zero
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0._wp
zfv_vw(ji,jj) = 0._wp
END_2D
ENDIF
!
DO jk = 1, jpk-2 ! divergence of advective fluxes
!
DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero
zfu_uw(ji,jj,jpk) = 0._wp ; zfv_vw(ji,jj,jpk) = 0._wp
zfu_uw(ji,jj, 1 ) = 0._wp ; zfv_vw(ji,jj, 1 ) = 0._wp
DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport at level k+1
zfw(ji,jj) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk+1)
END_2D
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
ENDIF
DO jk = 2, jpkm1 ! interior advective fluxes
DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport
zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk)
END_2D
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj ,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji ,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
END_2D
END DO
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! divergence of vertical momentum flux divergence
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) &
DO_2D( 0, 0, 0, 0 )
! ! vertical flux at level k+1
zzfu_kp1 = ( zfw(ji,jj) + zfw(ji+1,jj ) ) * ( puu(ji,jj,jk+1,Kmm) + puu(ji,jj,jk,Kmm) )
zzfv_kp1 = ( zfw(ji,jj) + zfw(ji ,jj+1) ) * ( pvv(ji,jj,jk+1,Kmm) + pvv(ji,jj,jk,Kmm) )
! ! divergence of vertical momentum flux
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj) - zzfu_kp1 ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_3D
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cen2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj) - zzfv_kp1 ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
! ! store vertical flux for next level calculation
zfu_uw(ji,jj) = zzfu_kp1
zfv_vw(ji,jj) = zzfv_kp1
END_2D
END DO
!
jk = jpkm1
DO_2D( 0, 0, 0, 0 )
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - zfu_uw(ji,jj) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - zfv_vw(ji,jj) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_zad, kt, Kmm )
DEALLOCATE( zu_trd, zv_trd )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cen2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
END SUBROUTINE dyn_adv_cen2
......
src/OCE/DYN/dynadv_ubs.F90
View file @ 4dbfc9bb
......@@ -6,6 +6,7 @@ MODULE dynadv_ubs
!!======================================================================
!! History : 2.0 ! 2006-08 (R. Benshila, L. Debreu) Original code
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -29,7 +30,8 @@ MODULE dynadv_ubs
REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp ! =1/4 quick ; =1/3 3rd order UBS
REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0 2nd order ; =1/32 4th order centred
PUBLIC dyn_adv_ubs ! routine called by step.F90
PUBLIC dyn_adv_ubs ! routine called by dynadv.F90
PUBLIC dyn_adv_ubs_hls1 ! routine called by dynadv.F90
!! * Substitutions
# include "do_loop_substitute.h90"
......@@ -41,7 +43,243 @@ MODULE dynadv_ubs
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_ubs ***
!!
!! ** Purpose : Compute the now momentum advection trend in flux form
!! and the general trend of the momentum equation.
!!
!! ** Method : The scheme is the one implemeted in ROMS. It depends
!! on two parameter gamma1 and gamma2. The former control the
!! upstream baised part of the scheme and the later the centred
!! part: gamma1 = 0 pure centered (no diffusive part)
!! = 1/4 Quick scheme
!! = 1/3 3rd order Upstream biased scheme
!! gamma2 = 0 2nd order finite differencing
!! = 1/32 4th order finite differencing
!! For stability reasons, the first term of the fluxes which cor-
!! responds to a second order centered scheme is evaluated using
!! the now velocity (centered in time) while the second term which
!! is the diffusive part of the scheme, is evaluated using the
!! before velocity (forward in time).
!! Default value (hard coded in the begining of the module) are
!! gamma1=1/3 and gamma2=1/32.
!!
!! In RK3 time stepping case, the optional arguments
!! (pau,pav,paw) are present. They are used as advective velocity
!! while the advected velocity remains (puu,pvv).
!!
!! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend
!!
!! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time-step and level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zzu, zui, zfuj, zl_u, zzfu_kp1 ! local scalars
REAL(wp) :: zzv, zvj, zfvi, zl_v, zzfv_kp1 ! - -
REAL(wp), DIMENSION(A2D(2)) :: zfu_t, zfu_f, zfu
REAL(wp), DIMENSION(A2D(2)) :: zfv_t, zfv_f, zfv
REAL(wp), DIMENSION(A2D(2),2) :: zlu_uu, zlu_uv
REAL(wp), DIMENSION(A2D(2),2) :: zlv_vv, zlv_vu
REAL(wp), DIMENSION(:,:,:) , POINTER :: zpt_u, zpt_v, zpt_w
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zu_trd, zv_trd
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
ALLOCATE( zu_trd(A2D(0),jpkm1), zv_trd(A2D(0),jpkm1) )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
zpt_u => pau(:,:,:)
zpt_v => pav(:,:,:)
zpt_w => paw(:,:,:)
ELSE ! MLF: advective velocity = (puu,pvv,ww)
zpt_u => puu(:,:,:,Kmm)
zpt_v => pvv(:,:,:,Kmm)
zpt_w => ww (:,:,: )
ENDIF
!
! ! =========================== !
DO jk = 1, jpkm1 ! Laplacian of the velocity !
! ! =========================== !
! ! horizontal volume fluxes
DO_2D( 2, 2, 2, 2 )
zfu(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
zfv(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
END_2D
!
DO_2D( 1, 1, 1, 1 ) ! laplacian
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility (north fold)
!! zlu_uu(ji,jj,1) = ( ( puu(ji+1,jj ,jk,Kbb) + puu(ji-1,jj ,jk,Kbb) ) - 2._wp * puu(ji,jj,jk,Kbb) ) * umask(ji,jj,jk)
!! zlv_vv(ji,jj,1) = ( ( pvv(ji ,jj+1,jk,Kbb) + pvv(ji ,jj-1,jk,Kbb) ) - 2._wp * pvv(ji,jj,jk,Kbb) ) * vmask(ji,jj,jk)
zlu_uu(ji,jj,1) = ( ( puu (ji+1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& + ( puu (ji-1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
& ) ) * umask(ji ,jj ,jk)
zlv_vv(ji,jj,1) = ( ( pvv (ji ,jj+1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& + ( pvv (ji ,jj-1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
& ) ) * vmask(ji ,jj ,jk)
zlu_uv(ji,jj,1) = ( puu(ji ,jj+1,jk,Kbb) - puu(ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
& - ( puu(ji ,jj ,jk,Kbb) - puu(ji ,jj-1,jk,Kbb) ) * fmask(ji ,jj-1,jk)
zlv_vu(ji,jj,1) = ( pvv(ji+1,jj ,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
& - ( pvv(ji ,jj ,jk,Kbb) - pvv(ji-1,jj ,jk,Kbb) ) * fmask(ji-1,jj ,jk)
!
!! zlu_uu(ji,jj,2) = ( ( zfu(ji+1,jj ) + zfu(ji-1,jj ) ) - 2._wp * zfu(ji,jj) ) * umask(ji ,jj ,jk)
!! zlv_vv(ji,jj,2) = ( ( zfv(ji ,jj+1) + zfv(ji ,jj-1) ) - 2._wp * zfv(ji,jj) ) * vmask(ji ,jj ,jk)
zlu_uu(ji,jj,2) = ( ( zfu(ji+1,jj ) - zfu(ji ,jj ) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& + ( zfu(ji-1,jj ) - zfu(ji ,jj ) &
& ) ) * umask(ji ,jj ,jk)
zlv_vv(ji,jj,2) = ( ( zfv(ji ,jj+1) - zfv(ji ,jj ) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& + ( zfv(ji ,jj-1) - zfv(ji ,jj ) &
& ) ) * vmask(ji ,jj ,jk)
zlu_uv(ji,jj,2) = ( zfu(ji ,jj+1) - zfu(ji ,jj ) ) * fmask(ji ,jj ,jk) &
& - ( zfu(ji ,jj ) - zfu(ji ,jj-1) ) * fmask(ji ,jj-1,jk)
zlv_vu(ji,jj,2) = ( zfv(ji+1,jj ) - zfv(ji ,jj ) ) * fmask(ji ,jj ,jk) &
& - ( zfv(ji ,jj ) - zfv(ji-1,jj ) ) * fmask(ji-1,jj ,jk)
END_2D
!
! ! ====================== !
! ! Horizontal advection !
! ! ====================== !
! ! horizontal volume fluxes
DO_2D( 1, 1, 1, 1 )
zfu(ji,jj) = 0.25_wp * zfu(ji,jj)
zfv(ji,jj) = 0.25_wp * zfv(ji,jj)
END_2D
!
DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes at T- and F-point
zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
!
IF( zui > 0 ) THEN ; zl_u = zlu_uu(ji ,jj,1)
ELSE ; zl_u = zlu_uu(ji+1,jj,1)
ENDIF
IF( zvj > 0 ) THEN ; zl_v = zlv_vv(ji,jj ,1)
ELSE ; zl_v = zlv_vv(ji,jj+1,1)
ENDIF
!
zfu_t(ji+1,jj ) = ( zfu(ji,jj) + zfu(ji+1,jj ) - gamma2 * ( zlu_uu(ji,jj,2) + zlu_uu(ji+1,jj ,2) ) ) &
& * ( zui - gamma1 * zl_u )
zfv_t(ji ,jj+1) = ( zfv(ji,jj) + zfv(ji ,jj+1) - gamma2 * ( zlv_vv(ji,jj,2) + zlv_vv(ji ,jj+1,2) ) ) &
& * ( zvj - gamma1 * zl_v )
!
zfuj = ( zfu(ji,jj) + zfu(ji ,jj+1) )
zfvi = ( zfv(ji,jj) + zfv(ji+1,jj ) )
IF( zfuj > 0 ) THEN ; zl_v = zlv_vu(ji ,jj,1)
ELSE ; zl_v = zlv_vu(ji+1,jj,1)
ENDIF
IF( zfvi > 0 ) THEN ; zl_u = zlu_uv(ji,jj ,1)
ELSE ; zl_u = zlu_uv(ji,jj+1,1)
ENDIF
!
zfv_f(ji ,jj ) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,2) + zlv_vu(ji+1,jj ,2) ) ) &
& * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) - gamma1 * zl_u )
zfu_f(ji ,jj ) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,2) + zlu_uv(ji ,jj+1,2) ) ) &
& * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) - gamma1 * zl_v )
END_2D
DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj) - zfu_t(ji,jj ) &
& + zfv_f(ji ,jj) - zfv_f(ji,jj-1) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ) - zfu_f(ji-1,jj) &
& + zfv_t(ji,jj+1) - zfv_t(ji ,jj) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
END DO
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
! ! ==================== !
! ! Vertical advection !
! ! ==================== !
!
#define zfu_uw zfu_t
#define zfv_vw zfv_t
#define zfw zfu
!
! ! surface vertical fluxes
!
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
ELSE ! non linear free: surface advective fluxes set to zero
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0._wp
zfv_vw(ji,jj) = 0._wp
END_2D
ENDIF
!
DO jk = 1, jpk-2 ! divergence of advective fluxes
!
DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport at level k+1
zfw(ji,jj) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk+1)
END_2D
DO_2D( 0, 0, 0, 0 )
! ! vertical flux at level k+1
zzfu_kp1 = ( zfw(ji,jj) + zfw(ji+1,jj ) ) * ( puu(ji,jj,jk+1,Kmm) + puu(ji,jj,jk,Kmm) )
zzfv_kp1 = ( zfw(ji,jj) + zfw(ji ,jj+1) ) * ( pvv(ji,jj,jk+1,Kmm) + pvv(ji,jj,jk,Kmm) )
! ! divergence of vertical momentum flux
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj) - zzfu_kp1 ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj) - zzfv_kp1 ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
! ! store vertical flux for next level calculation
zfu_uw(ji,jj) = zzfu_kp1
zfv_vw(ji,jj) = zzfv_kp1
END_2D
END DO
!
jk = jpkm1 ! compute last level (zzfu_kp1 = 0)
DO_2D( 0, 0, 0, 0 )
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - zfu_uw(ji,jj) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - zfv_vw(ji,jj) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_zad, kt, Kmm )
DEALLOCATE( zu_trd, zv_trd )
ENDIF
!
#undef zfu_uw
#undef zfv_vw
#undef zfw
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
END SUBROUTINE dyn_adv_ubs
SUBROUTINE dyn_adv_ubs_hls1( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_ubs ***
!!
......@@ -73,7 +311,6 @@ CONTAINS
!! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time-step and level indices
INTEGER , OPTIONAL , INTENT(in ) :: no_zad ! no vertical advection compotation
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
!
......@@ -89,7 +326,7 @@ CONTAINS
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
IF(lwp) WRITE(numout,*) 'dyn_adv_ubs_hls1 : UBS flux form momentum advection'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
......@@ -230,63 +467,44 @@ CONTAINS
! ! Vertical advection !
! ! ==================== !
!
! ! ======================== !
IF( PRESENT( no_zad ) ) THEN ! No vertical advection ! (except if linear free surface)
! ! ======================== ! ------
!
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,1,Kmm)
pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,1,Kmm)
END_2D
ENDIF
! ! =================== !
ELSE ! Vertical advection !
! ! =================== !
DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero
zfu_uw(ji,jj,jpk) = 0._wp
zfv_vw(ji,jj,jpk) = 0._wp
zfu_uw(ji,jj, 1 ) = 0._wp
zfv_vw(ji,jj, 1 ) = 0._wp
DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero
zfu_uw(ji,jj,jpk) = 0._wp
zfv_vw(ji,jj,jpk) = 0._wp
zfu_uw(ji,jj, 1 ) = 0._wp
zfv_vw(ji,jj, 1 ) = 0._wp
END_2D
IF( ln_linssh ) THEN ! constant volume : advection through the surface
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
IF( ln_linssh ) THEN ! constant volume : advection through the surface
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
ENDIF
DO jk = 2, jpkm1 ! interior fluxes
DO_2D( 0, 1, 0, 1 )
zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk)
END_2D
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
END_2D
END DO
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! divergence of vertical momentum flux divergence
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_3D
!
IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic
zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
ENDIF
DO jk = 2, jpkm1 ! interior fluxes
DO_2D( 0, 1, 0, 1 )
zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk)
END_2D
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
END_2D
END DO
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! divergence of vertical momentum flux divergence
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_3D
!
END SUBROUTINE dyn_adv_ubs
IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic
zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
END SUBROUTINE dyn_adv_ubs_hls1
!!==============================================================================
END MODULE dynadv_ubs
src/OCE/DYN/dynatf.F90
View file @ 4dbfc9bb
......@@ -328,29 +328,29 @@ CONTAINS
!
IF ( iom_use("utau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zutau(jpi,jpj))
ALLOCATE(zutau(A2D(0)))
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)
ELSE
CALL iom_put( "utau", utau(:,:) )
CALL iom_put( "utau", utau(A2D(0)) )
ENDIF
ENDIF
!
IF ( iom_use("vtau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zvtau(jpi,jpj))
ALLOCATE(zvtau(A2D(0)))
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)
ELSE
CALL iom_put( "vtau", vtau(:,:) )
CALL iom_put( "vtau", vtau(A2D(0)) )
ENDIF
ENDIF
!
......
src/OCE/DYN/dynatf_qco.F90
View file @ 4dbfc9bb
......@@ -245,29 +245,29 @@ CONTAINS
!
IF ( iom_use("utau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zutau(jpi,jpj))
ALLOCATE(zutau(A2D(0)))
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)
ELSE
CALL iom_put( "utau", utau(:,:) )
CALL iom_put( "utau", utau(A2D(0)) )
ENDIF
ENDIF
!
IF ( iom_use("vtau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zvtau(jpi,jpj))
ALLOCATE(zvtau(A2D(0)))
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)
ELSE
CALL iom_put( "vtau", vtau(:,:) )
CALL iom_put( "vtau", vtau(A2D(0)) )
ENDIF
ENDIF
!
......
src/OCE/DYN/dynkeg.F90
View file @ 4dbfc9bb
......@@ -6,7 +6,8 @@ MODULE dynkeg
!! History : 1.0 ! 1987-09 (P. Andrich, M.-A. Foujols) Original code
!! 7.0 ! 1997-05 (G. Madec) Split dynber into dynkeg and dynhpg
!! NEMO 1.0 ! 2002-07 (G. Madec) F90: Free form and module
!! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
!! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -27,7 +28,8 @@ MODULE dynkeg
IMPLICIT NONE
PRIVATE
PUBLIC dyn_keg ! routine called by step module
PUBLIC dyn_keg ! routine called by step module
PUBLIC dyn_keg_hls1 ! routine called by step module
INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme)
INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1 !: Hollingsworth et al., QJRMS, 1983
......@@ -44,6 +46,123 @@ MODULE dynkeg
CONTAINS
SUBROUTINE dyn_keg( kt, kscheme, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_keg ***
!!
!! ** Purpose : Compute the now momentum trend due to the horizontal
!! gradient of the horizontal kinetic energy and add it to the
!! general momentum trend.
!!
!! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that
!! conserve kinetic energy. Compute the now horizontal kinetic energy
!! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
!! * kscheme = nkeg_HW : Hollingsworth correction following
!! Arakawa (2001). The now horizontal kinetic energy is given by:
!! zhke = 1/6 [ mi-1( 2 * un^2 + ((u(j+1)+u(j-1))/2)^2 )
!! + mj-1( 2 * vn^2 + ((v(i+1)+v(i-1))/2)^2 ) ]
!!
!! Take its horizontal gradient and add it to the general momentum
!! trend.
!! u(rhs) = u(rhs) - 1/e1u di[ zhke ]
!! v(rhs) = v(rhs) - 1/e2v dj[ zhke ]
!!
!! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend
!! - send this trends to trd_dyn (l_trddyn=T) for post-processing
!!
!! ** References : Arakawa, A., International Geophysics 2001.
!! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: kscheme ! =0/1 type of KEG scheme
INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zu, zv ! local scalars
REAL(wp), DIMENSION(:,: ) , ALLOCATABLE :: zhke
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zu_trd, zv_trd
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_keg')
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
IF(lwp) WRITE(numout,*) '~~~~~~~'
ENDIF
ENDIF
!
IF( l_trddyn ) THEN ! Save the input trends
ALLOCATE( zu_trd(A2D(0),jpk), zv_trd(A2D(0),jpk) )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
SELECT CASE ( kscheme )
!
CASE ( nkeg_C2 ) !== Standard scheme ==!
ALLOCATE( zhke(A2D(1)) )
DO jk = 1, jpkm1
DO_2D( 0, 1, 0, 1 ) !* Horizontal kinetic energy at T-point
zu = puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) &
& + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm)
zv = pvv(ji ,jj-1,jk,Kmm) * pvv(ji ,jj-1,jk,Kmm) &
& + pvv(ji ,jj ,jk,Kmm) * pvv(ji ,jj ,jk,Kmm)
zhke(ji,jj) = 0.25_wp * ( zv + zu )
END_2D
!
DO_2D( 0, 0, 0, 0 ) !* grad( KE ) added to the general momentum trends
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj ) - zhke(ji,jj) ) * r1_e1u(ji,jj)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji ,jj+1) - zhke(ji,jj) ) * r1_e2v(ji,jj)
END_2D
END DO
DEALLOCATE( zhke )
!
CASE ( nkeg_HW ) !* Hollingsworth scheme
ALLOCATE( zhke(A2D(1)) )
DO jk = 1, jpkm1
DO_2D( 0, 1, 0, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zu = ( puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) &
& + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm) ) * 8._wp &
& + ( ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) * ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) &
& + ( puu(ji ,jj-1,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) * ( puu(ji ,jj-1,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) &
& ) ! bracket for halo 1 - halo 2 compatibility
zv = ( pvv(ji ,jj-1,jk,Kmm) * pvv(ji ,jj-1,jk,Kmm) &
& + pvv(ji ,jj ,jk,Kmm) * pvv(ji ,jj ,jk,Kmm) ) * 8._wp &
& + ( ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) * ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) &
& + ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) * ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) &
& ) ! bracket for halo 1 - halo 2 compatibility
zhke(ji,jj) = r1_48 * ( zv + zu )
END_2D
!
DO_2D( 0, 0, 0, 0 ) !* grad( KE ) added to the general momentum trends
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj ) - zhke(ji,jj) ) * r1_e1u(ji,jj)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji ,jj+1) - zhke(ji,jj) ) * r1_e2v(ji,jj)
END_2D
END DO
DEALLOCATE( zhke )
!
END SELECT
!
IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
DEALLOCATE( zu_trd, zv_trd )
ENDIF
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' keg - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
IF( ln_timing ) CALL timing_stop('dyn_keg')
!
END SUBROUTINE dyn_keg
SUBROUTINE dyn_keg_hls1( kt, kscheme, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_keg ***
!!
......@@ -86,7 +205,7 @@ CONTAINS
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
IF(lwp) WRITE(numout,*) 'dyn_keg_hls1 : kinetic energy gradient trend, scheme number=', kscheme
IF(lwp) WRITE(numout,*) '~~~~~~~'
ENDIF
ENDIF
......@@ -147,7 +266,7 @@ CONTAINS
!
IF( ln_timing ) CALL timing_stop('dyn_keg')
!
END SUBROUTINE dyn_keg
END SUBROUTINE dyn_keg_hls1
!!======================================================================
END MODULE dynkeg