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
Showing
with 2522 additions and 2317 deletions
This diff is collapsed.
......@@ -104,7 +104,7 @@ CONTAINS
!
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs
!! diagnostics
REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
REAL(wp), DIMENSION(A2D(0)) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
!!----------------------------------------------------------------------
!
IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
......@@ -133,8 +133,7 @@ CONTAINS
END DO
CALL icemax4D( ze_i , zei_max )
CALL icemax4D( ze_s , zes_max )
CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp, zes_max, 'T', 1._wp )
!
!
! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
......@@ -182,11 +181,11 @@ CONTAINS
DO jt = 1, icycle
! diagnostics
zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow
zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi
zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 )
zdiag_adv_mass(:,:) = SUM( pv_i (A2D(0),:) , dim=3 ) * rhoi + SUM( pv_s (A2D(0),:) , dim=3 ) * rhos &
& + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow
zdiag_adv_salt(:,:) = SUM( psv_i(A2D(0),:) , dim=3 ) * rhoi
zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(A2D(0),1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(A2D(0),1:nlay_s,:) , dim=4 ), dim=3 )
! record at_i before advection (for open water)
zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
......@@ -370,25 +369,25 @@ CONTAINS
ELSE
CALL lbc_lnk( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp )
ENDIF
CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp, pe_s, 'T', 1._wp )
!
!== Open water area ==!
zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
DO_2D( 0, 0, 0, 0 )
pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &
& - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
& - ( ( zudy(ji,jj) - zudy(ji-1,jj) ) & ! ad () for NP repro
& + ( zvdx(ji,jj) - zvdx(ji,jj-1) ) ) * r1_e1e2t(ji,jj) * zdt
END_2D
CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp )
!
! --- diagnostics --- !
diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow &
diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (A2D(0),:) , dim=3 ) * rhoi + SUM( pv_s (A2D(0),:) , dim=3 ) * rhos &
& + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow &
& - zdiag_adv_mass(:,:) ) * z1_dt
diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi &
diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(A2D(0),:) , dim=3 ) * rhoi &
& - zdiag_adv_salt(:,:) ) * z1_dt
diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) &
diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(A2D(0),1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(A2D(0),1:nlay_s,:) , dim=4 ), dim=3 ) &
& - zdiag_adv_heat(:,:) ) * z1_dt
!
! --- Ensure non-negative fields and in-bound thicknesses --- !
......@@ -518,7 +517,8 @@ CONTAINS
! thus we calculate the upstream solution and apply a limiter again
DO jl = 1, jpl
DO_2D( 0, 0, 0, 0 )
ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) )
ztra = - ( ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) & ! add () for NP repro
& + ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) )
!
zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
END_2D
......@@ -553,7 +553,8 @@ CONTAINS
! ---------------------------------
DO jl = 1, jpl
DO_2D( 0, 0, 0, 0 )
ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )
ztra = - ( ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) ) & ! add () for NP repro
& + ( zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) )
!
ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
END_2D
......@@ -645,10 +646,10 @@ CONTAINS
!
DO jl = 1, jpl !-- after tracer with upstream scheme
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) &
& + pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) &
& + ( pu (ji,jj ) - pu (ji-1,jj ) &
& + pv (ji,jj ) - pv (ji ,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk)
ztra = - ( ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) ) & ! add () for NP repro
& + ( pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) ) &
& + ( ( pu (ji,jj ) - pu (ji-1,jj ) ) &
& + ( pv (ji,jj ) - pv (ji ,jj-1 ) ) ) * pt(ji,jj,jl) * (1.-pamsk)
!
pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
END_2D
......@@ -912,7 +913,7 @@ CONTAINS
!
DO jl = 1, jpl
DO_2D( 1, 0, kloop, kloop )
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt(ji+1,jj,jl) + pt(ji,jj,jl) ) &
& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
END_2D
END DO
......@@ -922,7 +923,7 @@ CONTAINS
DO jl = 1, jpl
DO_2D( 1, 0, kloop, kloop )
zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt(ji+1,jj,jl) + pt(ji,jj,jl) ) &
& - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
END_2D
END DO
......@@ -934,9 +935,9 @@ CONTAINS
zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
zdx2 = e1u(ji,jj) * e1u(ji,jj)
!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) ) &
& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
& + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
& + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) ) &
& - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
END_2D
END DO
......@@ -948,9 +949,9 @@ CONTAINS
zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
zdx2 = e1u(ji,jj) * e1u(ji,jj)
!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) ) &
& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
& + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
& + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) ) &
& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
END_2D
END DO
......@@ -965,11 +966,11 @@ CONTAINS
zdx2 = e1u(ji,jj) * e1u(ji,jj)
!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
zdx4 = zdx2 * zdx2
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) ) &
& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
& + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
& + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) ) &
& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
& + r1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) &
& + r1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ((ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) ) &
& - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
END_2D
END DO
......@@ -983,7 +984,7 @@ CONTAINS
DO jl = 1, jpl
DO_2D( 1, 0, kloop, kloop )
IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt(ji+1,jj,jl) + pt(ji,jj,jl) ) &
& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
ENDIF
END_2D
......@@ -1050,7 +1051,7 @@ CONTAINS
CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
DO jl = 1, jpl
DO_2D( kloop, kloop, 1, 0 )
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) &
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) &
& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
END_2D
END DO
......@@ -1059,7 +1060,7 @@ CONTAINS
DO jl = 1, jpl
DO_2D( kloop, kloop, 1, 0 )
zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) &
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) &
& - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
END_2D
END DO
......@@ -1070,9 +1071,9 @@ CONTAINS
zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
zdy2 = e2v(ji,jj) * e2v(ji,jj)
!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) ) &
& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
& + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
& + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) ) &
& - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
END_2D
END DO
......@@ -1083,9 +1084,9 @@ CONTAINS
zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
zdy2 = e2v(ji,jj) * e2v(ji,jj)
!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) ) &
& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
& + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
& + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) ) &
& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
END_2D
END DO
......@@ -1100,11 +1101,11 @@ CONTAINS
zdy2 = e2v(ji,jj) * e2v(ji,jj)
!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
zdy4 = zdy2 * zdy2
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) ) &
& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
& + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
& + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) ) &
& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
& + r1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) &
& + r1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ((ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) ) &
& - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
END_2D
END DO
......@@ -1244,10 +1245,10 @@ CONTAINS
zneg = MAX( 0._wp, pfu_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj ,jl) ) &
& + MAX( 0._wp, pfv_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj-1,jl) )
!
zpos = zpos - (pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) &
& ) * ( 1. - pamsk )
zneg = zneg + (pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) &
& ) * ( 1. - pamsk )
zpos = zpos - ( pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) &
& + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk )
zneg = zneg + ( pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) &
& + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk )
!
! ! up & down beta terms
! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!)
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.