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 159 additions and 3315 deletions
tests/ICE_ADV2D/MY_SRC/usrdef_hgr.F90
View file @ 4dbfc9bb
......@@ -90,8 +90,8 @@ CONTAINS
#endif
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zti = REAL( mig0(ji), wp ) - 0.5_wp ! start at i=0.5 in the global grid without halos
ztj = REAL( mjg0(jj), wp ) - 0.5_wp ! start at j=0.5 in the global grid without halos
zti = REAL( mig(ji,0), wp ) - 0.5_wp ! start at i=0.5 in the global grid without halos
ztj = REAL( mjg(jj,0), wp ) - 0.5_wp ! start at j=0.5 in the global grid without halos
plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti
plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * ( zti + 0.5_wp )
......@@ -110,19 +110,19 @@ CONTAINS
!! clem: This can be used with a 1proc simulation but I think it breaks repro when >1procs are used
!! DO jj = 1, jpj
!! DO ji = 1, jpi
!! !!pe1t(ji,jj) = rn_dx * EXP( -0.8/REAL(jpiglo**2) * (mi0(ji)-REAL(jpiglo+1)*0.5)**2 ) ! gaussian shape
!! !!pe2t(ji,jj) = rn_dy * EXP( -0.8/REAL(jpjglo**2) * (mj0(jj)-REAL(jpjglo+1)*0.5)**2 ) ! gaussian shape
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) ) ! linear shape
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) ) ! linear shape
!! !!pe1t(ji,jj) = rn_dx * EXP( -0.8/REAL(jpiglo**2) * (mi0(ji,nn_hls)-REAL(jpiglo+1)*0.5)**2 ) ! gaussian shape
!! !!pe2t(ji,jj) = rn_dy * EXP( -0.8/REAL(jpjglo**2) * (mj0(jj,nn_hls)-REAL(jpjglo+1)*0.5)**2 ) ! gaussian shape
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji,nn_hls))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) ) ! linear shape
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj,nn_hls))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) ) ! linear shape
!! END DO
!! END DO
!!#if defined key_agrif
!! IF( .NOT. Agrif_Root() ) THEN ! only works if the zoom is positioned at the center of the parent grid
!! DO jj = 1, jpj
!! DO ji = 1, jpi
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) &
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji,nn_hls))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) &
!! & * REAL(jpiglo) / REAL(Agrif_Parent(jpiglo) * Agrif_Rhox()) ) ! factor to match parent grid
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) &
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj,nn_hls))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) &
!! & * REAL(jpjglo) / REAL(Agrif_Parent(jpjglo) * Agrif_Rhoy()) ) ! factor to match parent grid
!! END DO
!! END DO
......
tests/ICE_AGRIF/MY_SRC/usrdef_hgr.F90
View file @ 4dbfc9bb
......@@ -96,8 +96,8 @@ CONTAINS
ENDIF
#endif
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zti = REAL( mig0(ji)-1, wp ) ! start at i=0 in the global grid without halos
ztj = REAL( mjg0(jj)-1, wp ) ! start at j=0 in the global grid without halos
zti = REAL( mig(ji,0)-1, wp ) ! start at i=0 in the global grid without halos
ztj = REAL( mjg(jj,0)-1, wp ) ! start at j=0 in the global grid without halos
plamt(ji,jj) = roffsetx + rn_dx * 1.e-3 * ( zti - 0.5_wp )
plamu(ji,jj) = roffsetx + rn_dx * 1.e-3 * zti
......@@ -116,19 +116,19 @@ CONTAINS
!! clem: This can be used with a 1proc simulation but I think it breaks repro when >1procs are used
!! DO jj = 1, jpj
!! DO ji = 1, jpi
!! !!pe1t(ji,jj) = rn_dx * EXP( -0.8/REAL(jpiglo**2) * (mi0(ji)-REAL(jpiglo+1)*0.5)**2 ) ! gaussian shape
!! !!pe2t(ji,jj) = rn_dy * EXP( -0.8/REAL(jpjglo**2) * (mj0(jj)-REAL(jpjglo+1)*0.5)**2 ) ! gaussian shape
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) ) ! linear shape
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) ) ! linear shape
!! !!pe1t(ji,jj) = rn_dx * EXP( -0.8/REAL(jpiglo**2) * (mi0(ji,nn_hls)-REAL(jpiglo+1)*0.5)**2 ) ! gaussian shape
!! !!pe2t(ji,jj) = rn_dy * EXP( -0.8/REAL(jpjglo**2) * (mj0(jj,nn_hls)-REAL(jpjglo+1)*0.5)**2 ) ! gaussian shape
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji,nn_hls))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) ) ! linear shape
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj,nn_hls))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) ) ! linear shape
!! END DO
!! END DO
!!#if defined key_agrif
!! IF( .NOT. Agrif_Root() ) THEN ! only works if the zoom is positioned at the center of the parent grid
!! DO jj = 1, jpj
!! DO ji = 1, jpi
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) &
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji,nn_hls))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) &
!! & * REAL(jpiglo) / REAL(Agrif_Parent(jpiglo) * Agrif_Rhox()) ) ! factor to match parent grid
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) &
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj,nn_hls))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) &
!! & * REAL(jpjglo) / REAL(Agrif_Parent(jpjglo) * Agrif_Rhoy()) ) ! factor to match parent grid
!! END DO
!! END DO
......
tests/ICE_AGRIF/MY_SRC/usrdef_sbc.F90
View file @ 4dbfc9bb
......@@ -18,7 +18,7 @@ MODULE usrdef_sbc
USE sbc_ice ! Surface boundary condition: ice fields
USE phycst ! physical constants
USE ice, ONLY : at_i_b, a_i_b
USE icethd_dh ! for CALL ice_thd_snwblow
!! USE icethd_dh ! for CALL ice_thd_snwblow
USE sbc_phy, ONLY : pp_cldf
!
USE in_out_manager ! I/O manager
......@@ -33,6 +33,8 @@ MODULE usrdef_sbc
PUBLIC usrdef_sbc_ice_tau ! routine called by icestp.F90 for ice dynamics
PUBLIC usrdef_sbc_ice_flx ! routine called by icestp.F90 for ice thermo
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: usrdef_sbc.F90 14273 2021-01-06 10:57:45Z smasson $
......@@ -109,8 +111,8 @@ CONTAINS
!!
INTEGER :: jl
REAL(wp) :: zfr1, zfr2 ! local variables
REAL(wp), DIMENSION(jpi,jpj) :: zsnw ! snw distribution after wind blowing
REAL(wp), DIMENSION(jpi,jpj) :: ztri
REAL(wp), DIMENSION(A2D(0)) :: zsnw ! snw distribution after wind blowing
REAL(wp), DIMENSION(A2D(0)) :: ztri
!!---------------------------------------------------------------------
!
IF( kt==nit000 .AND. lwp) WRITE(numout,*)' usrdef_sbc_ice : ICE_AGRIF case: NO flux forcing'
......@@ -134,9 +136,9 @@ CONTAINS
emp_ice (:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw(:,:)
emp_oce (:,:) = emp_oce(:,:) - sprecip(:,:) * (1._wp - zsnw(:,:) )
qevap_ice(:,:,:) = 0._wp
qprec_ice(:,:) = rhos * ( sst_m(:,:) * rcpi - rLfus ) * tmask(:,:,1) ! in J/m3
qemp_oce (:,:) = - emp_oce(:,:) * sst_m(:,:) * rcp
qemp_ice (:,:) = sprecip(:,:) * zsnw * ( sst_m(:,:) * rcpi - rLfus ) * tmask(:,:,1) ! solid precip (only)
qprec_ice(:,:) = rhos * ( sst_m(A2D(0)) * rcpi - rLfus ) * smask0(:,:) ! in J/m3
qemp_oce (:,:) = - emp_oce(:,:) * sst_m(A2D(0)) * rcp
qemp_ice (:,:) = sprecip(:,:) * zsnw * ( sst_m(A2D(0)) * rcpi - rLfus ) * smask0(:,:) ! solid precip (only)
! total fluxes
emp_tot (:,:) = emp_ice + emp_oce
......@@ -148,11 +150,11 @@ CONTAINS
ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:) ! surface transmission when hi>10cm
!
DO jl = 1, jpl
WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp ) ! constant (ztri) when hi>10cm
WHERE ( phs(A2D(0),jl) <= 0._wp .AND. phi(A2D(0),jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(A2D(0),jl) * 10._wp ) )
ELSEWHERE( phs(A2D(0),jl) <= 0._wp .AND. phi(A2D(0),jl) >= 0.1_wp ) ! constant (ztri) when hi>10cm
qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:)
ELSEWHERE ! zero when hs>0
ELSEWHERE ! zero when hs>0
qtr_ice_top(:,:,jl) = 0._wp
END WHERE
ENDDO
......
tests/ICE_RHEO/EXPREF/file_def_nemo-oce.xml
View file @ 4dbfc9bb
......@@ -13,14 +13,17 @@
<file_group id="1ts" output_freq="1ts" output_level="10" enabled=".TRUE."/> <!-- 1 time step files -->
<file_group id="1h" output_freq="1h" output_level="10" enabled=".TRUE."> <!-- 1h files -->
<file id="file1" name_suffix="_grid_T" description="ocean T grid variables" >
<field field_ref="utau" name="sozotaux" />
<field field_ref="vtau" name="sometauy" />
</file>
<file id="file2" name_suffix="_grid_U" description="ocean U grid variables" >
<field field_ref="uoce" name="vozocrtx" />
<field field_ref="utau" name="sozotaux" />
</file>
<file id="file3" name_suffix="_grid_V" description="ocean V grid variables" >
<field field_ref="voce" name="vomecrty" />
<field field_ref="vtau" name="sometauy" />
</file>
</file_group>
......
tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90
View file @ 4dbfc9bb
......@@ -573,7 +573,7 @@ CONTAINS
& ) * zmsk00y(ji,jj)
ENDIF
!extra code for test case boundary conditions
IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
IF (mjg(jj,nn_hls)<25 .or. mjg(jj,nn_hls)>975 .or. mig(ji,nn_hls)<25 .or. mig(ji,nn_hls)>975) THEN
v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj))
END IF
......@@ -630,7 +630,7 @@ CONTAINS
& ) * zmsk00x(ji,jj)
ENDIF
!extra code for test case boundary conditions
IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
IF (mjg(jj,nn_hls)<25 .or. mjg(jj,nn_hls)>975 .or. mig(ji,nn_hls)<25 .or. mig(ji,nn_hls)>975) THEN
u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj))
END IF
......@@ -689,7 +689,7 @@ CONTAINS
& ) * zmsk00x(ji,jj)
ENDIF
!extra code for test case boundary conditions
IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
IF (mjg(jj,nn_hls)<25 .or. mjg(jj,nn_hls)>975 .or. mig(ji,nn_hls)<25 .or. mig(ji,nn_hls)>975) THEN
u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj))
END IF
END_2D
......@@ -745,7 +745,7 @@ CONTAINS
& ) * zmsk00y(ji,jj)
ENDIF
!extra code for test case boundary conditions
IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
IF (mjg(jj,nn_hls)<25 .or. mjg(jj,nn_hls)>975 .or. mig(ji,nn_hls)<25 .or. mig(ji,nn_hls)>975) THEN
v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj))
END IF
END_2D
......@@ -1048,7 +1048,7 @@ CONTAINS
zresm = 0._wp
DO_2D( 0, 0, 0, 0 )
! cut of the boundary of the box (forced velocities)
IF (mjg0(jj)>30 .AND. mjg0(jj)<=970 .AND. mig0(ji)>30 .AND. mig0(ji)<=970) THEN
IF (mjg(jj,0)>30 .AND. mjg(jj,0)<=970 .AND. mig(ji,0)>30 .AND. mig(ji,0)<=970) THEN
zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
& ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) )
ENDIF
......
tests/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90
View file @ 4dbfc9bb
......@@ -525,7 +525,7 @@ CONTAINS
& ) * zmsk00y(ji,jj)
ENDIF
!extra code for test case boundary conditions
IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
IF (mjg(jj,nn_hls)<25 .or. mjg(jj,nn_hls)>975 .or. mig(ji,nn_hls)<25 .or. mig(ji,nn_hls)>975) THEN
v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj))
END IF
END_2D
......@@ -581,7 +581,7 @@ CONTAINS
& ) * zmsk00x(ji,jj)
ENDIF
!extra code for test case boundary conditions
IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
IF (mjg(jj,nn_hls)<25 .or. mjg(jj,nn_hls)>975 .or. mig(ji,nn_hls)<25 .or. mig(ji,nn_hls)>975) THEN
u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj))
END IF
END_2D
......@@ -639,7 +639,7 @@ CONTAINS
& ) * zmsk00x(ji,jj)
ENDIF
!extra code for test case boundary conditions
IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
IF (mjg(jj,nn_hls)<25 .or. mjg(jj,nn_hls)>975 .or. mig(ji,nn_hls)<25 .or. mig(ji,nn_hls)>975) THEN
u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj))
END IF
END_2D
......@@ -695,7 +695,7 @@ CONTAINS
& ) * zmsk00y(ji,jj)
ENDIF
!extra code for test case boundary conditions
IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
IF (mjg(jj,nn_hls)<25 .or. mjg(jj,nn_hls)>975 .or. mig(ji,nn_hls)<25 .or. mig(ji,nn_hls)>975) THEN
v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj))
END IF
END_2D
......@@ -978,7 +978,7 @@ CONTAINS
zresm = 0._wp
DO_2D( 0, 0, 0, 0 )
! cut of the boundary of the box (forced velocities)
IF (mjg0(jj)>30 .AND. mjg0(jj)<=970 .AND. mig0(ji)>30 .AND. mig0(ji)<=970) THEN
IF (mjg(jj,0)>30 .AND. mjg(jj,0)<=970 .AND. mig(ji,0)>30 .AND. mig(ji,0)<=970) THEN
zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
& ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) )
ENDIF
......
tests/ICE_RHEO/MY_SRC/usrdef_hgr.F90
View file @ 4dbfc9bb
......@@ -98,17 +98,17 @@ CONTAINS
!! ==> EITHER 1) variable scale factors
!! clem: This can be used with a 1proc simulation but I think it breaks repro when >1procs are used
!! DO_2D( 1, 1, 1, 1 )
!! !!pe1t(ji,jj) = rn_dx * EXP( -0.8/REAL(jpiglo**2) * (mi0(ji)-REAL(jpiglo+1)*0.5)**2 ) ! gaussian shape
!! !!pe2t(ji,jj) = rn_dy * EXP( -0.8/REAL(jpjglo**2) * (mj0(jj)-REAL(jpjglo+1)*0.5)**2 ) ! gaussian shape
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) ) ! linear shape
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) ) ! linear shape
!! !!pe1t(ji,jj) = rn_dx * EXP( -0.8/REAL(jpiglo**2) * (mi0(ji,nn_hls)-REAL(jpiglo+1)*0.5)**2 ) ! gaussian shape
!! !!pe2t(ji,jj) = rn_dy * EXP( -0.8/REAL(jpjglo**2) * (mj0(jj,nn_hls)-REAL(jpjglo+1)*0.5)**2 ) ! gaussian shape
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji,nn_hls))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) ) ! linear shape
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj,nn_hls))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) ) ! linear shape
!! END_2D
!!#if defined key_agrif
!! IF( .NOT. Agrif_Root() ) THEN ! only works if the zoom is positioned at the center of the parent grid
!! DO_2D( 1, 1, 1, 1 )
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) &
!! pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji,nn_hls))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) &
!! & * REAL(jpiglo) / REAL(Agrif_Parent(jpiglo) * Agrif_Rhox()) ) ! factor to match parent grid
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) &
!! pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj,nn_hls))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) &
!! & * REAL(jpjglo) / REAL(Agrif_Parent(jpjglo) * Agrif_Rhoy()) ) ! factor to match parent grid
!! END_2D
!! ENDIF
......
tests/ICE_RHEO/MY_SRC/usrdef_nam.F90
View file @ 4dbfc9bb
......@@ -13,7 +13,6 @@ MODULE usrdef_nam
!! usr_def_nam : read user defined namelist and set global domain size
!! usr_def_hgr : initialize the horizontal mesh
!!----------------------------------------------------------------------
USE dom_oce , ONLY: nimpp , njmpp , Agrif_Root ! i- & j-indices of the local domain
USE par_oce ! ocean space and time domain
USE phycst ! physical constants
!
......
tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
View file @ 4dbfc9bb
......@@ -71,8 +71,8 @@ CONTAINS
!ij0 = 1 ; ij1 = 25 ! set boundary condition
!ii0 = 975 ; ii1 = 1000
!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)
! utau(ji,jj) = -utau_ice(ji,jj)
! vtau(ji,jj) = -vtau_ice(ji,jj)
! END DO
......@@ -108,7 +108,7 @@ CONTAINS
REAL(wp) :: zwndi_f , zwndj_f, zwnorm_f ! relative wind module and components at F-point
REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point
REAL(wp), DIMENSION(jpi,jpj) :: windu, windv ! wind components (idealised forcing)
REAL(wp), DIMENSION(A2D(0)) :: windu, windv ! wind components (idealised forcing)
REAL(wp), PARAMETER :: r_vfac = 1._wp ! relative velocity (make 0 for absolute velocity)
REAL(wp), PARAMETER :: Rwind = -0.8_wp ! ratio of wind components
REAL(wp), PARAMETER :: Umax = 15._wp ! maximum wind speed (m/s)
......@@ -122,40 +122,26 @@ CONTAINS
DO_2D( 0, 0, 0, 0 )
! wind spins up over 6 hours, factor 1000 to balance the units
windu(ji,jj) = Umax/sqrt(d*1000)*(d-2*mig(ji)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*min(kt*30./21600,1.)
windv(ji,jj) = Umax/sqrt(d*1000)*(d-2*mjg(jj)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*Rwind*min(kt*30./21600,1.)
windu(ji,jj) = Umax/SQRT(d*1000)*(d-2*mig(ji,nn_hls)*res) / &
& ((d-2*mig(ji,nn_hls)*res)**2+(d-2*mjg(jj,nn_hls)*res)**2*Rwind**2)**(1/4)*MIN(kt*30./21600,1.)
windv(ji,jj) = Umax/SQRT(d*1000)*(d-2*mjg(jj,nn_hls)*res) / &
& ((d-2*mig(ji,nn_hls)*res)**2+(d-2*mjg(jj,nn_hls)*res)**2*Rwind**2)**(1/4)*Rwind*MIN(kt*30./21600,1.)
END_2D
CALL lbc_lnk( 'usrdef_sbc', windu, 'U', -1., windv, 'V', -1. )
wndm_ice(:,:) = 0._wp !!gm brutal....
! ------------------------------------------------------------ !
! Wind module relative to the moving ice ( U10m - U_ice ) !
! Wind module and stress relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
! C-grid ice dynamics : U & V-points (same as ocean)
DO_2D( 0, 0, 0, 0 )
zwndi_t = ( windu(ji,jj) - r_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) )
zwndi_t = ( windu(ji,jj) - r_vfac * 0.5 * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) )
zwndj_t = ( windv(ji,jj) - r_vfac * 0.5 * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) )
!
wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
!
utau_ice(ji,jj) = zrhoa * Cd_atm * wndm_ice(ji,jj) * zwndi_t
vtau_ice(ji,jj) = zrhoa * Cd_atm * wndm_ice(ji,jj) * zwndj_t
END_2D
CALL lbc_lnk( 'usrdef_sbc', wndm_ice, 'T', 1. )
!!gm brutal....
utau_ice (:,:) = 0._wp
vtau_ice (:,:) = 0._wp
!!gm end
CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'T', -1., vtau_ice, 'T', -1., ldfull = .TRUE. )
! ------------------------------------------------------------ !
! Wind stress relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
! C-grid ice dynamics : U & V-points (same as ocean)
DO_2D( 0, 0, 0, 0 )
utau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) &
& * ( 0.5 * (windu(ji+1,jj) + windu(ji,jj) ) - r_vfac * u_ice(ji,jj) )
vtau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) &
& * ( 0.5 * (windv(ji,jj+1) + windv(ji,jj) ) - r_vfac * v_ice(ji,jj) )
END_2D
CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
!
END SUBROUTINE usrdef_sbc_ice_tau
......@@ -170,7 +156,7 @@ CONTAINS
REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness
!!
REAL(wp) :: zfr1, zfr2 ! local variables
REAL(wp), DIMENSION(jpi,jpj) :: zsnw ! snw distribution after wind blowing
REAL(wp), DIMENSION(A2D(0)) :: zsnw ! snw distribution after wind blowing
!!---------------------------------------------------------------------
!
IF( kt==nit000 .AND. lwp) WRITE(numout,*)' usrdef_sbc_ice : ICE_RHEO case: NO flux forcing'
......@@ -194,9 +180,9 @@ CONTAINS
emp_ice (:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw(:,:)
emp_oce (:,:) = emp_oce(:,:) - sprecip(:,:) * (1._wp - zsnw(:,:) )
qevap_ice(:,:,:) = 0._wp
qprec_ice(:,:) = rhos * ( sst_m(:,:) * rcpi - rLfus ) * tmask(:,:,1) ! in J/m3
qemp_oce (:,:) = - emp_oce(:,:) * sst_m(:,:) * rcp
qemp_ice (:,:) = sprecip(:,:) * zsnw * ( sst_m(:,:) * rcpi - rLfus ) * tmask(:,:,1) ! solid precip (only)
qprec_ice(:,:) = rhos * ( sst_m(A2D(0)) * rcpi - rLfus ) * smask0(:,:) ! in J/m3
qemp_oce (:,:) = - emp_oce(:,:) * sst_m(A2D(0)) * rcp
qemp_ice (:,:) = sprecip(:,:) * zsnw * ( sst_m(A2D(0)) * rcpi - rLfus ) * smask0(:,:) ! solid precip (only)
! total fluxes
emp_tot (:,:) = emp_ice + emp_oce
......@@ -207,11 +193,11 @@ CONTAINS
zfr1 = ( 0.18 * ( 1.0 - pp_cldf ) + 0.35 * pp_cldf ) ! transmission when hi>10cm
zfr2 = ( 0.82 * ( 1.0 - pp_cldf ) + 0.65 * pp_cldf ) ! zfr2 such that zfr1 + zfr2 to equal 1
!
WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm
WHERE ( phs(A2D(0),:) <= 0._wp .AND. phi(A2D(0),:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm
qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(A2D(0),:) * 10._wp ) )
ELSEWHERE( phs(A2D(0),:) <= 0._wp .AND. phi(A2D(0),:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm
qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
ELSEWHERE ! zero when hs>0
ELSEWHERE ! zero when hs>0
qtr_ice_top(:,:,:) = 0._wp
END WHERE
......
tests/ISOMIP+/EXPREF/namelist_cfg
View file @ 4dbfc9bb
......@@ -310,14 +310,19 @@ rn_Dt = 720.
!-----------------------------------------------------------------------
&nameos ! ocean Equation Of Seawater (default: NO selection)
!-----------------------------------------------------------------------
ln_leos = .true. ! = Use L-EOS (linear Eq.)
!
ln_seos = .true. ! = Use S-EOS (simplified Eq.)
! ! S-EOS coefficients (ln_seos=T):
! ! rd(T,S,Z)*rho0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS
! ! L-EOS coefficients (ln_seos=T):
! ! rd(T,S,Z)*rho0 = rho0*(-a0*dT+b0*dS)
! ! dT = T-rn_T0 ; dS = S-rn_S0
rn_T0 = -1. ! reference temperature
rn_S0 = 34.2 ! reference salinity
rn_a0 = 3.7330e-5 ! thermal expension coefficient
rn_b0 = 7.8430e-4 ! saline expension coefficient
rn_lambda1 = 0. ! cabbeling coeff in T^2 (=0 for linear eos)
rn_lambda2 = 0. ! cabbeling coeff in S^2 (=0 for linear eos)
rn_mu1 = 0. ! thermobaric coeff. in T (=0 for linear eos)
rn_mu2 = 0. ! thermobaric coeff. in S (=0 for linear eos)
rn_nu = 0. ! cabbeling coeff in T*S (=0 for linear eos)
/
!-----------------------------------------------------------------------
&namtra_adv ! advection scheme for tracer (default: NO selection)
......
tests/ISOMIP+/MY_SRC/dtatsd.F90
View file @ 4dbfc9bb
......@@ -33,8 +33,8 @@ MODULE dtatsd
LOGICAL , PUBLIC :: ln_tsd_init !: T & S data flag
LOGICAL , PUBLIC :: ln_tsd_dmp !: internal damping toward input data flag
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsd ! structure of input SST (file informations, fields read)
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsdini ! structure of input SST (file informations, fields read)
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsddmp ! structure of input SST (file informations, fields read)
!! * Substitutions
# include "do_loop_substitute.h90"
......@@ -121,16 +121,16 @@ CONTAINS
IF( ln_tsd_dmp ) THEN
!
ALLOCATE( sf_tsddmp(jpts), STAT=ierr0 )
ALLOCATE( sf_tsd(jpts), STAT=ierr0 )
IF( ierr0 > 0 ) THEN
CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsddmp structure' ) ; RETURN
CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' ) ; RETURN
ENDIF
!
! dmp file
ALLOCATE( sf_tsddmp(jp_tem)%fnow(jpi,jpj,jpk) , STAT=ierr0 )
IF( sn_dmpt%ln_tint ) ALLOCATE( sf_tsddmp(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )
ALLOCATE( sf_tsddmp(jp_sal)%fnow(jpi,jpj,jpk) , STAT=ierr2 )
IF( sn_dmps%ln_tint ) ALLOCATE( sf_tsddmp(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk) , STAT=ierr0 )
IF( sn_dmpt%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )
ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk) , STAT=ierr2 )
IF( sn_dmps%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
!
IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN
CALL ctl_stop( 'dta_tsd : unable to allocate T & S dmp data arrays' ) ; RETURN
......@@ -138,14 +138,14 @@ CONTAINS
!
! ! fill sf_tsd with sn_tem & sn_sal and control print
slf_i(jp_tem) = sn_dmpt ; slf_i(jp_sal) = sn_dmps
CALL fld_fill( sf_tsddmp, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity dmp data', 'namtsd', no_print )
CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity dmp data', 'namtsd', no_print )
!
ENDIF
!
END SUBROUTINE dta_tsd_init
SUBROUTINE dta_tsd( kt, cddta, ptsd )
SUBROUTINE dta_tsd( kt, ptsd, cddta )
!!----------------------------------------------------------------------
!! *** ROUTINE dta_tsd ***
!!
......@@ -159,45 +159,43 @@ CONTAINS
!!
!! ** Action : ptsd T-S data on medl mesh and interpolated at time-step kt
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step
CHARACTER(LEN=3) , INTENT(in ) :: cddta ! dmp or ini
REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts), INTENT( out) :: ptsd ! T & S data
INTEGER , INTENT(in ) :: kt ! ocean time-step
REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts), INTENT( out) :: ptsd ! T & S data
CHARACTER(len=*), OPTIONAL , INTENT(in ) :: cddta ! force the initialization when tradmp is used
!
INTEGER :: ji, jj, jk, jl, jkk ! dummy loop indicies
INTEGER :: ik, il0, il1, ii0, ii1, ij0, ij1 ! local integers
REAL(wp):: zl, zi ! local scalars
LOGICAL :: ll_tsdini
REAL(wp), DIMENSION(jpk) :: ztp, zsp ! 1D workspace
!!----------------------------------------------------------------------
!
ll_tsdini = .FALSE.
IF( PRESENT(cddta) ) ll_tsdini = .TRUE.
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only for the full domain
IF( ln_tile ) CALL dom_tile_stop( ldhold=.TRUE. ) ! Use full domain
SELECT CASE(cddta)
CASE('ini')
IF( ll_tsdini ) THEN
CALL fld_read( kt, 1, sf_tsdini ) !== read T & S data at kt time step ==!
CASE('dmp')
CALL fld_read( kt, 1, sf_tsddmp ) !== read T & S data at kt time step ==!
CASE DEFAULT
CALL ctl_stop('STOP', 'dta_tsd: cddta case unknown')
END SELECT
ELSE
CALL fld_read( kt, 1, sf_tsd ) !== read T & S data at kt time step ==!
ENDIF
IF( ln_tile ) CALL dom_tile_start( ldhold=.TRUE. ) ! Revert to tile domain
ENDIF
!
SELECT CASE(cddta)
CASE('ini')
IF( ll_tsdini ) THEN
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
ptsd(ji,jj,jk,jp_tem) = sf_tsdini(jp_tem)%fnow(ji,jj,jk) ! NO mask
ptsd(ji,jj,jk,jp_sal) = sf_tsdini(jp_sal)%fnow(ji,jj,jk)
END_3D
CASE('dmp')
ELSE
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
ptsd(ji,jj,jk,jp_tem) = sf_tsddmp(jp_tem)%fnow(ji,jj,jk) ! NO mask
ptsd(ji,jj,jk,jp_sal) = sf_tsddmp(jp_sal)%fnow(ji,jj,jk)
ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk) ! NO mask
ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk)
END_3D
CASE DEFAULT
CALL ctl_stop('STOP', 'dta_tsd: cddta case unknown')
END SELECT
ENDIF
!
IF( ln_sco ) THEN !== s- or mixed s-zps-coordinate ==!
!
......@@ -261,8 +259,7 @@ CONTAINS
!
ENDIF
!
SELECT CASE(cddta)
CASE('ini')
IF( ll_tsdini ) THEN
! !== deallocate T & S structure ==!
! (data used only for initialisation)
IF(lwp) WRITE(numout,*) 'dta_tsd: deallocte T & S arrays as they are only use to initialize the run'
......@@ -272,7 +269,7 @@ CONTAINS
IF( sf_tsdini(jp_sal)%ln_tint ) DEALLOCATE( sf_tsdini(jp_sal)%fdta )
DEALLOCATE( sf_tsdini ) ! the structure itself
!
END SELECT
ENDIF
!
END SUBROUTINE dta_tsd
......
tests/ISOMIP+/MY_SRC/eosbn2.F90 deleted 100644 → 0
View file @ 0b9e2937
MODULE eosbn2
!!==============================================================================
!! *** MODULE eosbn2 ***
!! Equation Of Seawater : in situ density - Brunt-Vaisala frequency
!!==============================================================================
!! History : OPA ! 1989-03 (O. Marti) Original code
!! 6.0 ! 1994-07 (G. Madec, M. Imbard) add bn2
!! 6.0 ! 1994-08 (G. Madec) Add Jackett & McDougall eos
!! 7.0 ! 1996-01 (G. Madec) statement function for e3
!! 8.1 ! 1997-07 (G. Madec) density instead of volumic mass
!! - ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure gradient
!! 8.2 ! 2001-09 (M. Ben Jelloul) bugfix on linear eos
!! NEMO 1.0 ! 2002-10 (G. Madec) add eos_init
!! - ! 2002-11 (G. Madec, A. Bozec) partial step, eos_insitu_2d
!! - ! 2003-08 (G. Madec) F90, free form
!! 3.0 ! 2006-08 (G. Madec) add tfreez function (now eos_fzp function)
!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA
!! - ! 2010-10 (G. Nurser, G. Madec) add alpha/beta used in ldfslp
!! 3.7 ! 2012-03 (F. Roquet, G. Madec) add primitive of alpha and beta used in PE computation
!! - ! 2012-05 (F. Roquet) add Vallis and original JM95 equation of state
!! - ! 2013-04 (F. Roquet, G. Madec) add eos_rab, change bn2 computation and reorganize the module
!! - ! 2014-09 (F. Roquet) add TEOS-10, S-EOS, and modify EOS-80
!! - ! 2015-06 (P.A. Bouttier) eos_fzp functions changed to subroutines for AGRIF
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! eos : generic interface of the equation of state
!! eos_insitu : Compute the in situ density
!! eos_insitu_pot: Compute the insitu and surface referenced potential volumic mass
!! eos_insitu_2d : Compute the in situ density for 2d fields
!! bn2 : compute the Brunt-Vaisala frequency
!! eos_pt_from_ct: compute the potential temperature from the Conservative Temperature
!! eos_rab : generic interface of in situ thermal/haline expansion ratio
!! eos_rab_3d : compute in situ thermal/haline expansion ratio
!! eos_rab_2d : compute in situ thermal/haline expansion ratio for 2d fields
!! eos_fzp_2d : freezing temperature for 2d fields
!! eos_fzp_0d : freezing temperature for scalar
!! eos_init : set eos parameters (namelist)
!!----------------------------------------------------------------------
USE dom_oce ! ocean space and time domain
USE domutl, ONLY : is_tile
USE phycst ! physical constants
USE stopar ! Stochastic T/S fluctuations
USE stopts ! Stochastic T/S fluctuations
!
USE in_out_manager ! I/O manager
USE lib_mpp ! MPP library
USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
USE prtctl ! Print control
USE lbclnk ! ocean lateral boundary conditions
USE timing ! Timing
IMPLICIT NONE
PRIVATE
! !! * Interface
INTERFACE eos
MODULE PROCEDURE eos_insitu_New, eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d
END INTERFACE
!
INTERFACE eos_rab
MODULE PROCEDURE rab_3d, rab_2d, rab_0d
END INTERFACE
!
INTERFACE eos_fzp
MODULE PROCEDURE eos_fzp_2d, eos_fzp_0d
END INTERFACE
!
PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules
PUBLIC bn2 ! called by step module
PUBLIC eos_rab ! called by ldfslp, zdfddm, trabbl
PUBLIC eos_pt_from_ct ! called by sbcssm
PUBLIC eos_fzp ! called by traadv_cen2 and sbcice_... modules
PUBLIC eos_pen ! used for pe diagnostics in trdpen module
PUBLIC eos_init ! called by istate module
! !!** Namelist nameos **
LOGICAL , PUBLIC :: ln_TEOS10
LOGICAL , PUBLIC :: ln_EOS80
LOGICAL , PUBLIC :: ln_SEOS
LOGICAL , PUBLIC :: ln_LEOS ! determine if linear eos is used
! Parameters
LOGICAL , PUBLIC :: l_useCT ! =T in ln_TEOS10=T (i.e. use eos_pt_from_ct to compute sst_m), =F otherwise
INTEGER , PUBLIC :: neos ! Identifier for equation of state used
INTEGER , PARAMETER :: np_teos10 = -1 ! parameter for using TEOS10
INTEGER , PARAMETER :: np_eos80 = 0 ! parameter for using EOS80
INTEGER , PARAMETER :: np_seos = 1 ! parameter for using Simplified Equation of state
INTEGER , PARAMETER :: np_leos = 2 ! parameter for using linear equation of state (ISOMIP+)
! !!! simplified eos coefficients (default value: Vallis 2006)
REAL(wp), PUBLIC :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff.
REAL(wp), PUBLIC :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff.
REAL(wp) :: rn_lambda1 = 5.9520e-2_wp ! cabbeling coeff. in T^2
REAL(wp) :: rn_lambda2 = 5.4914e-4_wp ! cabbeling coeff. in S^2
REAL(wp) :: rn_mu1 = 1.4970e-4_wp ! thermobaric coeff. in T
REAL(wp) :: rn_mu2 = 1.1090e-5_wp ! thermobaric coeff. in S
REAL(wp) :: rn_nu = 2.4341e-3_wp ! cabbeling coeff. in theta*salt
! TEOS10/EOS80 parameters
REAL(wp) :: r1_S0, r1_T0, r1_Z0, rdeltaS
! EOS parameters
REAL(wp) :: EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600
REAL(wp) :: EOS010 , EOS110 , EOS210 , EOS310 , EOS410 , EOS510
REAL(wp) :: EOS020 , EOS120 , EOS220 , EOS320 , EOS420
REAL(wp) :: EOS030 , EOS130 , EOS230 , EOS330
REAL(wp) :: EOS040 , EOS140 , EOS240
REAL(wp) :: EOS050 , EOS150
REAL(wp) :: EOS060
REAL(wp) :: EOS001 , EOS101 , EOS201 , EOS301 , EOS401
REAL(wp) :: EOS011 , EOS111 , EOS211 , EOS311
REAL(wp) :: EOS021 , EOS121 , EOS221
REAL(wp) :: EOS031 , EOS131
REAL(wp) :: EOS041
REAL(wp) :: EOS002 , EOS102 , EOS202
REAL(wp) :: EOS012 , EOS112
REAL(wp) :: EOS022
REAL(wp) :: EOS003 , EOS103
REAL(wp) :: EOS013
! ALPHA parameters
REAL(wp) :: ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500
REAL(wp) :: ALP010 , ALP110 , ALP210 , ALP310 , ALP410
REAL(wp) :: ALP020 , ALP120 , ALP220 , ALP320
REAL(wp) :: ALP030 , ALP130 , ALP230
REAL(wp) :: ALP040 , ALP140
REAL(wp) :: ALP050
REAL(wp) :: ALP001 , ALP101 , ALP201 , ALP301
REAL(wp) :: ALP011 , ALP111 , ALP211
REAL(wp) :: ALP021 , ALP121
REAL(wp) :: ALP031
REAL(wp) :: ALP002 , ALP102
REAL(wp) :: ALP012
REAL(wp) :: ALP003
! BETA parameters
REAL(wp) :: BET000 , BET100 , BET200 , BET300 , BET400 , BET500
REAL(wp) :: BET010 , BET110 , BET210 , BET310 , BET410
REAL(wp) :: BET020 , BET120 , BET220 , BET320
REAL(wp) :: BET030 , BET130 , BET230
REAL(wp) :: BET040 , BET140
REAL(wp) :: BET050
REAL(wp) :: BET001 , BET101 , BET201 , BET301
REAL(wp) :: BET011 , BET111 , BET211
REAL(wp) :: BET021 , BET121
REAL(wp) :: BET031
REAL(wp) :: BET002 , BET102
REAL(wp) :: BET012
REAL(wp) :: BET003
! PEN parameters
REAL(wp) :: PEN000 , PEN100 , PEN200 , PEN300 , PEN400
REAL(wp) :: PEN010 , PEN110 , PEN210 , PEN310
REAL(wp) :: PEN020 , PEN120 , PEN220
REAL(wp) :: PEN030 , PEN130
REAL(wp) :: PEN040
REAL(wp) :: PEN001 , PEN101 , PEN201
REAL(wp) :: PEN011 , PEN111
REAL(wp) :: PEN021
REAL(wp) :: PEN002 , PEN102
REAL(wp) :: PEN012
! ALPHA_PEN parameters
REAL(wp) :: APE000 , APE100 , APE200 , APE300
REAL(wp) :: APE010 , APE110 , APE210
REAL(wp) :: APE020 , APE120
REAL(wp) :: APE030
REAL(wp) :: APE001 , APE101
REAL(wp) :: APE011
REAL(wp) :: APE002
! BETA_PEN parameters
REAL(wp) :: BPE000 , BPE100 , BPE200 , BPE300
REAL(wp) :: BPE010 , BPE110 , BPE210
REAL(wp) :: BPE020 , BPE120
REAL(wp) :: BPE030
REAL(wp) :: BPE001 , BPE101
REAL(wp) :: BPE011
REAL(wp) :: BPE002
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: eosbn2.F90 10425 2018-12-19 21:54:16Z smasson $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE eos_insitu_New( pts, Knn, prd )
!!----------------------------------------------------------------------
!! *** ROUTINE eos_insitu_New ***
!!
!! ** Purpose : Compute the in situ density (ratio rho/rho0) from
!! potential temperature and salinity using an equation of state
!! selected in the nameos namelist
!!
!! ** Method : prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0
!! with prd in situ density anomaly no units
!! t TEOS10: CT or EOS80: PT Celsius
!! s TEOS10: SA or EOS80: SP TEOS10: g/kg or EOS80: psu
!! z depth meters
!! rho in situ density kg/m^3
!! rho0 reference density kg/m^3
!!
!! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
!! Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg
!!
!! ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z).
!! Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu
!!
!! ln_seos : simplified equation of state
!! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0
!! linear case function of T only: rn_alpha<>0, other coefficients = 0
!! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0
!! Vallis like equation: use default values of coefficients
!!
!! ln_leos : linear ISOMIP equation of state
!! prd(t,s,z) = ( -a0*(T-T0) + b0*(S-S0) ) / rho0
!! setup for ISOMIP linear eos
!!
!! ** Action : compute prd , the in situ density (no units)
!!
!! References : Roquet et al, Ocean Modelling, in preparation (2014)
!! Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
!! TEOS-10 Manual, 2010
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:,:,:,:), INTENT(in ) :: pts ! T-S
INTEGER , INTENT(in ) :: Knn ! time-level
REAL(wp), DIMENSION(:,:,: ), INTENT( out) :: prd ! in situ density
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos-insitu')
!
SELECT CASE( neos )
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
DO_3D(nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
!
zh = gdept(ji,jj,jk,Knn) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem,Knn) * r1_T0 ! temperature
zs = SQRT( ABS( pts(ji,jj,jk,jp_sal,Knn) + rdeltaS ) * r1_S0 ) ! square root salinity
ztm = tmask(ji,jj,jk) ! tmask
!
zn3 = EOS013*zt &
& + EOS103*zs+EOS003
!
zn2 = (EOS022*zt &
& + EOS112*zs+EOS012)*zt &
& + (EOS202*zs+EOS102)*zs+EOS002
!
zn1 = (((EOS041*zt &
& + EOS131*zs+EOS031)*zt &
& + (EOS221*zs+EOS121)*zs+EOS021)*zt &
& + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &
& + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
!
zn0 = (((((EOS060*zt &
& + EOS150*zs+EOS050)*zt &
& + (EOS240*zs+EOS140)*zs+EOS040)*zt &
& + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
& + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
& + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
& + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
!
zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
!
prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked)
!
END_3D
!
CASE( np_seos ) !== simplified EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem,Knn) - 10._wp
zs = pts (ji,jj,jk,jp_sal,Knn) - 35._wp
zh = gdept(ji,jj,jk,Knn)
ztm = tmask(ji,jj,jk)
!
zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt &
& + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs &
& - rn_nu * zt * zs
!
prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked)
END_3D
!
CASE( np_leos ) !== linear ISOMIP EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem,Knn) - (-1._wp)
zs = pts (ji,jj,jk,jp_sal,Knn) - 34.2_wp
zh = gdept(ji,jj,jk, Knn)
ztm = tmask(ji,jj,jk)
!
zn = rho0 * ( - rn_a0 * zt + rn_b0 * zs )
!
prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked)
END_3D
!
END SELECT
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu : ', kdim=jpk )
!
IF( ln_timing ) CALL timing_stop('eos-insitu')
!
END SUBROUTINE eos_insitu_New
SUBROUTINE eos_insitu( pts, prd, pdep )
!!
REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
! ! 2 : salinity [psu]
REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-]
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m]
!!
CALL eos_insitu_t( pts, is_tile(pts), prd, is_tile(prd), pdep, is_tile(pdep) )
END SUBROUTINE eos_insitu
SUBROUTINE eos_insitu_t( pts, ktts, prd, ktrd, pdep, ktdep )
!!----------------------------------------------------------------------
!! *** ROUTINE eos_insitu ***
!!
!! ** Purpose : Compute the in situ density (ratio rho/rho0) from
!! potential temperature and salinity using an equation of state
!! selected in the nameos namelist
!!
!! ** Method : prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0
!! with prd in situ density anomaly no units
!! t TEOS10: CT or EOS80: PT Celsius
!! s TEOS10: SA or EOS80: SP TEOS10: g/kg or EOS80: psu
!! z depth meters
!! rho in situ density kg/m^3
!! rho0 reference density kg/m^3
!!
!! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
!! Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg
!!
!! ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z).
!! Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu
!!
!! ln_seos : simplified equation of state
!! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0
!! linear case function of T only: rn_alpha<>0, other coefficients = 0
!! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0
!! Vallis like equation: use default values of coefficients
!!
!! ln_leos : linear ISOMIP equation of state
!! prd(t,s,z) = ( -a0*(T-T0) + b0*(S-S0) ) / rho0
!! setup for ISOMIP linear eos
!!
!! ** Action : compute prd , the in situ density (no units)
!!
!! References : Roquet et al, Ocean Modelling, in preparation (2014)
!! Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
!! TEOS-10 Manual, 2010
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ktts, ktrd, ktdep
REAL(wp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
! ! 2 : salinity [psu]
REAL(wp), DIMENSION(A2D_T(ktrd) ,JPK ), INTENT( out) :: prd ! in situ density [-]
REAL(wp), DIMENSION(A2D_T(ktdep),JPK ), INTENT(in ) :: pdep ! depth [m]
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos-insitu')
!
SELECT CASE( neos )
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
!
zh = pdep(ji,jj,jk) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
ztm = tmask(ji,jj,jk) ! tmask
!
zn3 = EOS013*zt &
& + EOS103*zs+EOS003
!
zn2 = (EOS022*zt &
& + EOS112*zs+EOS012)*zt &
& + (EOS202*zs+EOS102)*zs+EOS002
!
zn1 = (((EOS041*zt &
& + EOS131*zs+EOS031)*zt &
& + (EOS221*zs+EOS121)*zs+EOS021)*zt &
& + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &
& + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
!
zn0 = (((((EOS060*zt &
& + EOS150*zs+EOS050)*zt &
& + (EOS240*zs+EOS140)*zs+EOS040)*zt &
& + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
& + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
& + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
& + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
!
zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
!
prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked)
!
END_3D
!
CASE( np_seos ) !== simplified EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - 10._wp
zs = pts (ji,jj,jk,jp_sal) - 35._wp
zh = pdep (ji,jj,jk)
ztm = tmask(ji,jj,jk)
!
zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt &
& + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs &
& - rn_nu * zt * zs
!
prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked)
END_3D
!
CASE( np_leos ) !== linear ISOMIP EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - (-1._wp)
zs = pts (ji,jj,jk,jp_sal) - 34.2_wp
zh = pdep (ji,jj,jk)
ztm = tmask(ji,jj,jk)
!
zn = rho0 * ( - rn_a0 * zt + rn_b0 * zs )
!
prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked)
END_3D
!
END SELECT
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu : ', kdim=jpk )
!
IF( ln_timing ) CALL timing_stop('eos-insitu')
!
END SUBROUTINE eos_insitu_t
SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
!!
REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
! ! 2 : salinity [psu]
REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-]
REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prhop ! potential density (surface referenced)
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m]
!!
CALL eos_insitu_pot_t( pts, is_tile(pts), prd, is_tile(prd), prhop, is_tile(prhop), pdep, is_tile(pdep) )
END SUBROUTINE eos_insitu_pot
SUBROUTINE eos_insitu_pot_t( pts, ktts, prd, ktrd, prhop, ktrhop, pdep, ktdep )
!!----------------------------------------------------------------------
!! *** ROUTINE eos_insitu_pot ***
!!
!! ** Purpose : Compute the in situ density (ratio rho/rho0) and the
!! potential volumic mass (Kg/m3) from potential temperature and
!! salinity fields using an equation of state selected in the
!! namelist.
!!
!! ** Action : - prd , the in situ density (no units)
!! - prhop, the potential volumic mass (Kg/m3)
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ktts, ktrd, ktrhop, ktdep
REAL(wp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
! ! 2 : salinity [psu]
REAL(wp), DIMENSION(A2D_T(ktrd) ,JPK ), INTENT( out) :: prd ! in situ density [-]
REAL(wp), DIMENSION(A2D_T(ktrhop),JPK ), INTENT( out) :: prhop ! potential density (surface referenced)
REAL(wp), DIMENSION(A2D_T(ktdep) ,JPK ), INTENT(in ) :: pdep ! depth [m]
!
INTEGER :: ji, jj, jk, jsmp ! dummy loop indices
INTEGER :: jdof
REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos-pot')
!
SELECT CASE ( neos )
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
! Stochastic equation of state
IF ( ln_sto_eos ) THEN
ALLOCATE(zn0_sto(1:2*nn_sto_eos))
ALLOCATE(zn_sto(1:2*nn_sto_eos))
ALLOCATE(zsign(1:2*nn_sto_eos))
DO jsmp = 1, 2*nn_sto_eos, 2
zsign(jsmp) = 1._wp
zsign(jsmp+1) = -1._wp
END DO
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
!
! compute density (2*nn_sto_eos) times:
! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
! (2) for t-dt, s-ds (with the opposite fluctuation)
DO jsmp = 1, nn_sto_eos*2
jdof = (jsmp + 1) / 2
zh = pdep(ji,jj,jk) * r1_Z0 ! depth
zt = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0 ! temperature
zstemp = pts (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)
zs = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 ) ! square root salinity
ztm = tmask(ji,jj,jk) ! tmask
!
zn3 = EOS013*zt &
& + EOS103*zs+EOS003
!
zn2 = (EOS022*zt &
& + EOS112*zs+EOS012)*zt &
& + (EOS202*zs+EOS102)*zs+EOS002
!
zn1 = (((EOS041*zt &
& + EOS131*zs+EOS031)*zt &
& + (EOS221*zs+EOS121)*zs+EOS021)*zt &
& + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &
& + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
!
zn0_sto(jsmp) = (((((EOS060*zt &
& + EOS150*zs+EOS050)*zt &
& + (EOS240*zs+EOS140)*zs+EOS040)*zt &
& + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
& + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
& + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
& + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
!
zn_sto(jsmp) = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)
END DO
!
! compute stochastic density as the mean of the (2*nn_sto_eos) densities
prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp
DO jsmp = 1, nn_sto_eos*2
prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface
!
prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rho0 - 1._wp ) ! density anomaly (masked)
END DO
prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
prd (ji,jj,jk) = 0.5_wp * prd (ji,jj,jk) * ztm / nn_sto_eos
END_3D
DEALLOCATE(zn0_sto,zn_sto,zsign)
! Non-stochastic equation of state
ELSE
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
!
zh = pdep(ji,jj,jk) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
ztm = tmask(ji,jj,jk) ! tmask
!
zn3 = EOS013*zt &
& + EOS103*zs+EOS003
!
zn2 = (EOS022*zt &
& + EOS112*zs+EOS012)*zt &
& + (EOS202*zs+EOS102)*zs+EOS002
!
zn1 = (((EOS041*zt &
& + EOS131*zs+EOS031)*zt &
& + (EOS221*zs+EOS121)*zs+EOS021)*zt &
& + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &
& + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
!
zn0 = (((((EOS060*zt &
& + EOS150*zs+EOS050)*zt &
& + (EOS240*zs+EOS140)*zs+EOS040)*zt &
& + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
& + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
& + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
& + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
!
zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
!
prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface
!
prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked)
END_3D
ENDIF
CASE( np_seos ) !== simplified EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - 10._wp
zs = pts (ji,jj,jk,jp_sal) - 35._wp
zh = pdep (ji,jj,jk)
ztm = tmask(ji,jj,jk)
! ! potential density referenced at the surface
zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt &
& + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs &
& - rn_nu * zt * zs
prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
! ! density anomaly (masked)
zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
prd(ji,jj,jk) = zn * r1_rho0 * ztm
!
END_3D
!
CASE( np_leos ) !== linear ISOMIP EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - (-1._wp)
zs = pts (ji,jj,jk,jp_sal) - 34.2_wp
zh = pdep (ji,jj,jk)
ztm = tmask(ji,jj,jk)
! ! potential density referenced at the surface
zn = rho0 * ( - rn_a0 * zt + rn_b0 * zs )
prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
! ! density anomaly (masked)
prd(ji,jj,jk) = zn * r1_rho0 * ztm
!
END_3D
!
END SELECT
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', &
& tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk )
!
IF( ln_timing ) CALL timing_stop('eos-pot')
!
END SUBROUTINE eos_insitu_pot_t
SUBROUTINE eos_insitu_2d( pts, pdep, prd )
!!
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
! ! 2 : salinity [psu]
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pdep ! depth [m]
REAL(wp), DIMENSION(:,:) , INTENT( out) :: prd ! in situ density
!!
CALL eos_insitu_2d_t( pts, is_tile(pts), pdep, is_tile(pdep), prd, is_tile(prd) )
END SUBROUTINE eos_insitu_2d
SUBROUTINE eos_insitu_2d_t( pts, ktts, pdep, ktdep, prd, ktrd )
!!----------------------------------------------------------------------
!! *** ROUTINE eos_insitu_2d ***
!!
!! ** Purpose : Compute the in situ density (ratio rho/rho0) from
!! potential temperature and salinity using an equation of state
!! selected in the nameos namelist. * 2D field case
!!
!! ** Action : - prd , the in situ density (no units) (unmasked)
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ktts, ktdep, ktrd
REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
! ! 2 : salinity [psu]
REAL(wp), DIMENSION(A2D_T(ktdep) ), INTENT(in ) :: pdep ! depth [m]
REAL(wp), DIMENSION(A2D_T(ktrd) ), INTENT( out) :: prd ! in situ density
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos2d')
!
prd(:,:) = 0._wp
!
SELECT CASE( neos )
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
zh = pdep(ji,jj) * r1_Z0 ! depth
zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature
zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
!
zn3 = EOS013*zt &
& + EOS103*zs+EOS003
!
zn2 = (EOS022*zt &
& + EOS112*zs+EOS012)*zt &
& + (EOS202*zs+EOS102)*zs+EOS002
!
zn1 = (((EOS041*zt &
& + EOS131*zs+EOS031)*zt &
& + (EOS221*zs+EOS121)*zs+EOS021)*zt &
& + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &
& + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
!
zn0 = (((((EOS060*zt &
& + EOS150*zs+EOS050)*zt &
& + (EOS240*zs+EOS140)*zs+EOS040)*zt &
& + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
& + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
& + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
& + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
!
zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
!
prd(ji,jj) = zn * r1_rho0 - 1._wp ! unmasked in situ density anomaly
!
END_2D
!
CASE( np_seos ) !== simplified EOS ==!
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
zt = pts (ji,jj,jp_tem) - 10._wp
zs = pts (ji,jj,jp_sal) - 35._wp
zh = pdep (ji,jj) ! depth at the partial step level
!
zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt &
& + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs &
& - rn_nu * zt * zs
!
prd(ji,jj) = zn * r1_rho0 ! unmasked in situ density anomaly
!
END_2D
!
CASE( np_leos ) !== ISOMIP EOS ==!
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
zt = pts (ji,jj,jp_tem) - (-1._wp)
zs = pts (ji,jj,jp_sal) - 34.2_wp
zh = pdep (ji,jj) ! depth at the partial step level
!
zn = rho0 * ( - rn_a0 * zt + rn_b0 * zs )
!
prd(ji,jj) = zn * r1_rho0 ! unmasked in situ density anomaly
!
END_2D
!
!
END SELECT
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
!
IF( ln_timing ) CALL timing_stop('eos2d')
!
END SUBROUTINE eos_insitu_2d_t
SUBROUTINE eos_insitu_pot_2d( pts, prhop )
!!
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
! ! 2 : salinity [psu]
REAL(wp), DIMENSION(:,:) , INTENT( out) :: prhop ! potential density (surface referenced)
!!
CALL eos_insitu_pot_2d_t( pts, is_tile(pts), prhop, is_tile(prhop) )
END SUBROUTINE eos_insitu_pot_2d
SUBROUTINE eos_insitu_pot_2d_t( pts, ktts, prhop, ktrhop )
!!----------------------------------------------------------------------
!! *** ROUTINE eos_insitu_pot ***
!!
!! ** Purpose : Compute the in situ density (ratio rho/rho0) and the
!! potential volumic mass (Kg/m3) from potential temperature and
!! salinity fields using an equation of state selected in the
!! namelist.
!!
!! ** Action :
!! - prhop, the potential volumic mass (Kg/m3)
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ktts, ktrhop
REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
! ! 2 : salinity [psu]
REAL(wp), DIMENSION(A2D_T(ktrhop) ), INTENT( out) :: prhop ! potential density (surface referenced)
!
INTEGER :: ji, jj, jk, jsmp ! dummy loop indices
INTEGER :: jdof
REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos-pot')
!
SELECT CASE ( neos )
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature
zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
ztm = tmask(ji,jj,1) ! tmask
!
zn0 = (((((EOS060*zt &
& + EOS150*zs+EOS050)*zt &
& + (EOS240*zs+EOS140)*zs+EOS040)*zt &
& + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
& + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
& + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
& + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
!
!
prhop(ji,jj) = zn0 * ztm ! potential density referenced at the surface
!
END_2D
CASE( np_seos ) !== simplified EOS ==!
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zt = pts (ji,jj,jp_tem) - 10._wp
zs = pts (ji,jj,jp_sal) - 35._wp
ztm = tmask(ji,jj,1)
! ! potential density referenced at the surface
zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt &
& + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs &
& - rn_nu * zt * zs
prhop(ji,jj) = ( rho0 + zn ) * ztm
!
END_2D
!
CASE( np_leos ) !== ISOMIP EOS ==!
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
zt = pts (ji,jj,jp_tem) - (-1._wp)
zs = pts (ji,jj,jp_sal) - 34.2_wp
!zh = pdep (ji,jj) ! depth at the partial step level
!
zn = rho0 * ( - rn_a0 * zt + rn_b0 * zs )
!
prhop(ji,jj) = zn * r1_rho0 ! unmasked in situ density anomaly
!
END_2D
!
END SELECT
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=prhop, clinfo1=' pot: ', kdim=1 )
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=prhop, clinfo1=' eos-pot: ' )
!
IF( ln_timing ) CALL timing_stop('eos-pot')
!
END SUBROUTINE eos_insitu_pot_2d_t
SUBROUTINE rab_3d( pts, pab, Kmm )
!!
INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity
REAL(wp), DIMENSION(:,:,:,:), INTENT( out) :: pab ! thermal/haline expansion ratio
!!
CALL rab_3d_t( pts, is_tile(pts), pab, is_tile(pab), Kmm )
END SUBROUTINE rab_3d
SUBROUTINE rab_3d_t( pts, ktts, pab, ktab, Kmm )
!!----------------------------------------------------------------------
!! *** ROUTINE rab_3d ***
!!
!! ** Purpose : Calculates thermal/haline expansion ratio at T-points
!!
!! ** Method : calculates alpha / beta at T-points
!!
!! ** Action : - pab : thermal/haline expansion ratio at T-points
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: Kmm ! time level index
INTEGER , INTENT(in ) :: ktts, ktab
REAL(wp), DIMENSION(A2D_T(ktts),JPK,JPTS), INTENT(in ) :: pts ! pot. temperature & salinity
REAL(wp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT( out) :: pab ! thermal/haline expansion ratio
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('rab_3d')
!
SELECT CASE ( neos )
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
!
zh = gdept(ji,jj,jk,Kmm) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
ztm = tmask(ji,jj,jk) ! tmask
!
! alpha
zn3 = ALP003
!
zn2 = ALP012*zt + ALP102*zs+ALP002
!
zn1 = ((ALP031*zt &
& + ALP121*zs+ALP021)*zt &
& + (ALP211*zs+ALP111)*zs+ALP011)*zt &
& + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
!
zn0 = ((((ALP050*zt &
& + ALP140*zs+ALP040)*zt &
& + (ALP230*zs+ALP130)*zs+ALP030)*zt &
& + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt &
& + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt &
& + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
!
zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
!
pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm
!
! beta
zn3 = BET003
!
zn2 = BET012*zt + BET102*zs+BET002
!
zn1 = ((BET031*zt &
& + BET121*zs+BET021)*zt &
& + (BET211*zs+BET111)*zs+BET011)*zt &
& + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
!
zn0 = ((((BET050*zt &
& + BET140*zs+BET040)*zt &
& + (BET230*zs+BET130)*zs+BET030)*zt &
& + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt &
& + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt &
& + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
!
zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
!
pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm
!
END_3D
!
CASE( np_seos ) !== simplified EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0)
zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point
ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask
!
zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm ! alpha
!
zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm ! beta
!
END_3D
!
CASE( np_leos ) !== linear ISOMIP EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - (-1._wp)
zs = pts (ji,jj,jk,jp_sal) - 34.2_wp ! abs. salinity anomaly (s-S0)
zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point
ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask
!
zn = rn_a0 * rho0
pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm ! alpha
!
zn = rn_b0 * rho0
pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm ! beta
!
END_3D
!
CASE DEFAULT
WRITE(ctmp1,*) ' bad flag value for neos = ', neos
CALL ctl_stop( 'rab_3d:', ctmp1 )
!
END SELECT
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', &
& tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk )
!
IF( ln_timing ) CALL timing_stop('rab_3d')
!
END SUBROUTINE rab_3d_t
SUBROUTINE rab_2d( pts, pdep, pab, Kmm )
!!
INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pdep ! depth [m]
REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pab ! thermal/haline expansion ratio
!!
CALL rab_2d_t(pts, is_tile(pts), pdep, is_tile(pdep), pab, is_tile(pab), Kmm)
END SUBROUTINE rab_2d
SUBROUTINE rab_2d_t( pts, ktts, pdep, ktdep, pab, ktab, Kmm )
!!----------------------------------------------------------------------
!! *** ROUTINE rab_2d ***
!!
!! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked)
!!
!! ** Action : - pab : thermal/haline expansion ratio at T-points
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: Kmm ! time level index
INTEGER , INTENT(in ) :: ktts, ktdep, ktab
REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! pot. temperature & salinity
REAL(wp), DIMENSION(A2D_T(ktdep) ), INTENT(in ) :: pdep ! depth [m]
REAL(wp), DIMENSION(A2D_T(ktab),JPTS), INTENT( out) :: pab ! thermal/haline expansion ratio
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('rab_2d')
!
pab(:,:,:) = 0._wp
!
SELECT CASE ( neos )
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
zh = pdep(ji,jj) * r1_Z0 ! depth
zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature
zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
!
! alpha
zn3 = ALP003
!
zn2 = ALP012*zt + ALP102*zs+ALP002
!
zn1 = ((ALP031*zt &
& + ALP121*zs+ALP021)*zt &
& + (ALP211*zs+ALP111)*zs+ALP011)*zt &
& + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
!
zn0 = ((((ALP050*zt &
& + ALP140*zs+ALP040)*zt &
& + (ALP230*zs+ALP130)*zs+ALP030)*zt &
& + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt &
& + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt &
& + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
!
zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
!
pab(ji,jj,jp_tem) = zn * r1_rho0
!
! beta
zn3 = BET003
!
zn2 = BET012*zt + BET102*zs+BET002
!
zn1 = ((BET031*zt &
& + BET121*zs+BET021)*zt &
& + (BET211*zs+BET111)*zs+BET011)*zt &
& + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
!
zn0 = ((((BET050*zt &
& + BET140*zs+BET040)*zt &
& + (BET230*zs+BET130)*zs+BET030)*zt &
& + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt &
& + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt &
& + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
!
zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
!
pab(ji,jj,jp_sal) = zn / zs * r1_rho0
!
!
END_2D
!
CASE( np_seos ) !== simplified EOS ==!
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
zt = pts (ji,jj,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0)
zs = pts (ji,jj,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
zh = pdep (ji,jj) ! depth at the partial step level
!
zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
pab(ji,jj,jp_tem) = zn * r1_rho0 ! alpha
!
zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
pab(ji,jj,jp_sal) = zn * r1_rho0 ! beta
!
END_2D
!
CASE( np_leos ) !== linear ISOMIP EOS ==!
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
zt = pts (ji,jj,jp_tem) - (-1._wp) ! pot. temperature anomaly (t-T0)
zs = pts (ji,jj,jp_sal) - 34.2_wp ! abs. salinity anomaly (s-S0)
zh = pdep (ji,jj) ! depth at the partial step level
!
zn = rn_a0 * rho0
pab(ji,jj,jp_tem) = zn * r1_rho0 ! alpha
!
zn = rn_b0 * rho0
pab(ji,jj,jp_sal) = zn * r1_rho0 ! beta
!
END_2D
!
CASE DEFAULT
WRITE(ctmp1,*) ' bad flag value for neos = ', neos
CALL ctl_stop( 'rab_2d:', ctmp1 )
!
END SELECT
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', &
& tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )
!
IF( ln_timing ) CALL timing_stop('rab_2d')
!
END SUBROUTINE rab_2d_t
SUBROUTINE rab_0d( pts, pdep, pab, Kmm )
!!----------------------------------------------------------------------
!! *** ROUTINE rab_0d ***
!!
!! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked)
!!
!! ** Action : - pab : thermal/haline expansion ratio at T-points
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(jpts) , INTENT(in ) :: pts ! pot. temperature & salinity
REAL(wp), INTENT(in ) :: pdep ! depth [m]
REAL(wp), DIMENSION(jpts) , INTENT( out) :: pab ! thermal/haline expansion ratio
!
REAL(wp) :: zt , zh , zs ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('rab_0d')
!
pab(:) = 0._wp
!
SELECT CASE ( neos )
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
!
zh = pdep * r1_Z0 ! depth
zt = pts (jp_tem) * r1_T0 ! temperature
zs = SQRT( ABS( pts(jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
!
! alpha
zn3 = ALP003
!
zn2 = ALP012*zt + ALP102*zs+ALP002
!
zn1 = ((ALP031*zt &
& + ALP121*zs+ALP021)*zt &
& + (ALP211*zs+ALP111)*zs+ALP011)*zt &
& + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
!
zn0 = ((((ALP050*zt &
& + ALP140*zs+ALP040)*zt &
& + (ALP230*zs+ALP130)*zs+ALP030)*zt &
& + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt &
& + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt &
& + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
!
zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
!
pab(jp_tem) = zn * r1_rho0
!
! beta
zn3 = BET003
!
zn2 = BET012*zt + BET102*zs+BET002
!
zn1 = ((BET031*zt &
& + BET121*zs+BET021)*zt &
& + (BET211*zs+BET111)*zs+BET011)*zt &
& + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
!
zn0 = ((((BET050*zt &
& + BET140*zs+BET040)*zt &
& + (BET230*zs+BET130)*zs+BET030)*zt &
& + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt &
& + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt &
& + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
!
zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
!
pab(jp_sal) = zn / zs * r1_rho0
!
!
!
CASE( np_seos ) !== simplified EOS ==!
!
zt = pts(jp_tem) - 10._wp ! pot. temperature anomaly (t-T0)
zs = pts(jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
zh = pdep ! depth at the partial step level
!
zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
pab(jp_tem) = zn * r1_rho0 ! alpha
!
zn = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
pab(jp_sal) = zn * r1_rho0 ! beta
!
CASE( np_leos ) !== linear ISOMIP EOS ==!
!
zt = pts(jp_tem) - (-1._wp) ! pot. temperature anomaly (t-T0)
zs = pts(jp_sal) - 34.2_wp ! abs. salinity anomaly (s-S0)
zh = pdep ! depth at the partial step level
!
zn = rn_a0 * rho0
pab(jp_tem) = zn * r1_rho0 ! alpha
!
zn = rn_b0 * rho0
pab(jp_sal) = zn * r1_rho0 ! beta
!
CASE DEFAULT
WRITE(ctmp1,*) ' bad flag value for neos = ', neos
CALL ctl_stop( 'rab_0d:', ctmp1 )
!
END SELECT
!
IF( ln_timing ) CALL timing_stop('rab_0d')
!
END SUBROUTINE rab_0d
SUBROUTINE bn2( pts, pab, pn2, Kmm )
!!
INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu]
REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pab ! thermal/haline expansion coef. [Celsius-1,psu-1]
REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: pn2 ! Brunt-Vaisala frequency squared [1/s^2]
!!
CALL bn2_t( pts, pab, is_tile(pab), pn2, is_tile(pn2), Kmm )
END SUBROUTINE bn2
SUBROUTINE bn2_t( pts, pab, ktab, pn2, ktn2, Kmm )
!!----------------------------------------------------------------------
!! *** ROUTINE bn2 ***
!!
!! ** Purpose : Compute the local Brunt-Vaisala frequency at the
!! time-step of the input arguments
!!
!! ** Method : pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w
!! where alpha and beta are given in pab, and computed on T-points.
!! N.B. N^2 is set one for all to zero at jk=1 in istate module.
!!
!! ** Action : pn2 : square of the brunt-vaisala frequency at w-point
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: Kmm ! time level index
INTEGER , INTENT(in ) :: ktab, ktn2
REAL(wp), DIMENSION(jpi,jpj, jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu]
REAL(wp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT(in ) :: pab ! thermal/haline expansion coef. [Celsius-1,psu-1]
REAL(wp), DIMENSION(A2D_T(ktn2),JPK ), INTENT( out) :: pn2 ! Brunt-Vaisala frequency squared [1/s^2]
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zaw, zbw, zrw ! local scalars
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('bn2')
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 ) ! interior points only (2=< jk =< jpkm1 ); surface and bottom value set to zero one for all in istate.F90
zrw = ( gdepw(ji,jj,jk ,Kmm) - gdept(ji,jj,jk,Kmm) ) &
& / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) )
!
zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw
zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
!
pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) &
& - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) &
& / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
END_3D
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', kdim=jpk )
!
IF( ln_timing ) CALL timing_stop('bn2')
!
END SUBROUTINE bn2_t
FUNCTION eos_pt_from_ct( ctmp, psal ) RESULT( ptmp )
!!----------------------------------------------------------------------
!! *** ROUTINE eos_pt_from_ct ***
!!
!! ** Purpose : Compute pot.temp. from cons. temp. [Celsius]
!!
!! ** Method : rational approximation (5/3th order) of TEOS-10 algorithm
!! checkvalue: pt=20.02391895 Celsius for sa=35.7g/kg, ct=20degC
!!
!! Reference : TEOS-10, UNESCO
!! Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC)
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ctmp ! Cons. Temp [Celsius]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]
! Leave result array automatic rather than making explicitly allocated
REAL(wp), DIMENSION(jpi,jpj) :: ptmp ! potential temperature [Celsius]
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zt , zs , ztm ! local scalars
REAL(wp) :: zn , zd ! local scalars
REAL(wp) :: zdeltaS , z1_S0 , z1_T0
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos_pt_from_ct')
!
zdeltaS = 5._wp
z1_S0 = 0.875_wp/35.16504_wp
z1_T0 = 1._wp/40._wp
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
zt = ctmp (ji,jj) * z1_T0
zs = SQRT( ABS( psal(ji,jj) + zdeltaS ) * z1_S0 )
ztm = tmask(ji,jj,1)
!
zn = ((((-2.1385727895e-01_wp*zt &
& - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt &
& + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt &
& + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt &
& + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs &
& +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt &
& + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs &
& -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
!
zd = (2.0035003456_wp*zt &
& -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt &
& + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
!
ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
!
END_2D
!
IF( ln_timing ) CALL timing_stop('eos_pt_from_ct')
!
END FUNCTION eos_pt_from_ct
SUBROUTINE eos_fzp_2d( psal, ptf, pdep )
!!
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ), OPTIONAL :: pdep ! depth [m]
REAL(wp), DIMENSION(:,:) , INTENT(out ) :: ptf ! freezing temperature [Celsius]
!!
CALL eos_fzp_2d_t( psal, ptf, is_tile(ptf), pdep )
END SUBROUTINE eos_fzp_2d
SUBROUTINE eos_fzp_2d_t( psal, ptf, kttf, pdep )
!!----------------------------------------------------------------------
!! *** ROUTINE eos_fzp ***
!!
!! ** Purpose : Compute the freezing point temperature [Celsius]
!!
!! ** Method : UNESCO freezing point (ptf) in Celsius is given by
!! ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
!! checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
!!
!! Reference : UNESCO tech. papers in the marine science no. 28. 1978
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kttf
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: psal ! salinity [psu]
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ), OPTIONAL :: pdep ! depth [m]
REAL(wp), DIMENSION(A2D_T(kttf)), INTENT(out ) :: ptf ! freezing temperature [Celsius]
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zt, zs, z1_S0 ! local scalars
!!----------------------------------------------------------------------
!
SELECT CASE ( neos )
!
CASE ( np_teos10, np_seos ) !== CT,SA (TEOS-10 and S-EOS formulations) ==!
!
z1_S0 = 1._wp / 35.16504_wp
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 ) ! square root salinity
ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
& - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
END_2D
ptf(:,:) = ptf(:,:) * psal(:,:)
!
IF( PRESENT( pdep ) ) ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
!
CASE ( np_eos80 ) !== PT,SP (UNESCO formulation) ==!
!
ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) ) &
& - 2.154996e-4_wp * psal(:,:) ) * psal(:,:)
!
IF( PRESENT( pdep ) ) ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
!
CASE DEFAULT
WRITE(ctmp1,*) ' bad flag value for neos = ', neos
CALL ctl_stop( 'eos_fzp_2d:', ctmp1 )
!
END SELECT
!
END SUBROUTINE eos_fzp_2d_t
SUBROUTINE eos_fzp_0d( psal, ptf, pdep )
!!----------------------------------------------------------------------
!! *** ROUTINE eos_fzp ***
!!
!! ** Purpose : Compute the freezing point temperature [Celsius]
!!
!! ** Method : UNESCO freezing point (ptf) in Celsius is given by
!! ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z
!! checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
!!
!! Reference : UNESCO tech. papers in the marine science no. 28. 1978
!!----------------------------------------------------------------------
REAL(wp), INTENT(in ) :: psal ! salinity [psu]
REAL(wp), INTENT(in ), OPTIONAL :: pdep ! depth [m]
REAL(wp), INTENT(out) :: ptf ! freezing temperature [Celsius]
!
REAL(wp) :: zs ! local scalars
!!----------------------------------------------------------------------
!
SELECT CASE ( neos )
!
CASE ( np_teos10, np_seos ) !== CT,SA (TEOS-10 and S-EOS formulations) ==!
!
zs = SQRT( ABS( psal ) / 35.16504_wp ) ! square root salinity
ptf = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
& - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
ptf = ptf * psal
!
IF( PRESENT( pdep ) ) ptf = ptf - 7.53e-4 * pdep
!
CASE ( np_eos80 ) !== PT,SP (UNESCO formulation) ==!
!
ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal ) &
& - 2.154996e-4_wp * psal ) * psal
!
IF( PRESENT( pdep ) ) ptf = ptf - 7.53e-4 * pdep
!
CASE DEFAULT
WRITE(ctmp1,*) ' bad flag value for neos = ', neos
CALL ctl_stop( 'eos_fzp_0d:', ctmp1 )
!
END SELECT
!
END SUBROUTINE eos_fzp_0d
SUBROUTINE eos_pen( pts, pab_pe, ppen, Kmm )
!!----------------------------------------------------------------------
!! *** ROUTINE eos_pen ***
!!
!! ** Purpose : Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points
!!
!! ** Method : PE is defined analytically as the vertical
!! primitive of EOS times -g integrated between 0 and z>0.
!! pen is the nonlinear bsq-PE anomaly: pen = ( PE - rho0 gz ) / rho0 gz - rd
!! = 1/z * /int_0^z rd dz - rd
!! where rd is the density anomaly (see eos_rhd function)
!! ab_pe are partial derivatives of PE anomaly with respect to T and S:
!! ab_pe(1) = - 1/(rho0 gz) * dPE/dT + drd/dT = - d(pen)/dT
!! ab_pe(2) = 1/(rho0 gz) * dPE/dS + drd/dS = d(pen)/dS
!!
!! ** Action : - pen : PE anomaly given at T-points
!! : - pab_pe : given at T-points
!! pab_pe(:,:,:,jp_tem) is alpha_pe
!! pab_pe(:,:,:,jp_sal) is beta_pe
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab_pe ! alpha_pe and beta_pe
REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: ppen ! potential energy anomaly
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs , ztm ! local scalars
REAL(wp) :: zn , zn0, zn1, zn2 ! - -
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('eos_pen')
!
SELECT CASE ( neos )
!
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
!
zh = gdept(ji,jj,jk,Kmm) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
ztm = tmask(ji,jj,jk) ! tmask
!
! potential energy non-linear anomaly
zn2 = (PEN012)*zt &
& + PEN102*zs+PEN002
!
zn1 = ((PEN021)*zt &
& + PEN111*zs+PEN011)*zt &
& + (PEN201*zs+PEN101)*zs+PEN001
!
zn0 = ((((PEN040)*zt &
& + PEN130*zs+PEN030)*zt &
& + (PEN220*zs+PEN120)*zs+PEN020)*zt &
& + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt &
& + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000
!
zn = ( zn2 * zh + zn1 ) * zh + zn0
!
ppen(ji,jj,jk) = zn * zh * r1_rho0 * ztm
!
! alphaPE non-linear anomaly
zn2 = APE002
!
zn1 = (APE011)*zt &
& + APE101*zs+APE001
!
zn0 = (((APE030)*zt &
& + APE120*zs+APE020)*zt &
& + (APE210*zs+APE110)*zs+APE010)*zt &
& + ((APE300*zs+APE200)*zs+APE100)*zs+APE000
!
zn = ( zn2 * zh + zn1 ) * zh + zn0
!
pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm
!
! betaPE non-linear anomaly
zn2 = BPE002
!
zn1 = (BPE011)*zt &
& + BPE101*zs+BPE001
!
zn0 = (((BPE030)*zt &
& + BPE120*zs+BPE020)*zt &
& + (BPE210*zs+BPE110)*zs+BPE010)*zt &
& + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000
!
zn = ( zn2 * zh + zn1 ) * zh + zn0
!
pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm
!
END_3D
!
CASE( np_seos ) !== Vallis (2006) simplified EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts(ji,jj,jk,jp_tem) - 10._wp ! temperature anomaly (t-T0)
zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point
ztm = tmask(ji,jj,jk) ! tmask
zn = 0.5_wp * zh * r1_rho0 * ztm
! ! Potential Energy
ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
! ! alphaPE
pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn
pab_pe(ji,jj,jk,jp_sal) = rn_b0 * rn_mu2 * zn
!
END_3D
!
CASE( np_leos ) !== linear ISOMIP EOS ==!
!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts(ji,jj,jk,jp_tem) - (-1._wp) ! temperature anomaly (t-T0)
zs = pts (ji,jj,jk,jp_sal) - 34.2_wp ! abs. salinity anomaly (s-S0)
zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point
ztm = tmask(ji,jj,jk) ! tmask
zn = 0.5_wp * zh * r1_rho0 * ztm
! ! Potential Energy
ppen(ji,jj,jk) = 0.
! ! alphaPE
pab_pe(ji,jj,jk,jp_tem) = 0.
pab_pe(ji,jj,jk,jp_sal) = 0.
!
END_3D
!
CASE DEFAULT
WRITE(ctmp1,*) ' bad flag value for neos = ', neos
CALL ctl_stop( 'eos_pen:', ctmp1 )
!
END SELECT
!
IF( ln_timing ) CALL timing_stop('eos_pen')
!
END SUBROUTINE eos_pen
SUBROUTINE eos_init
!!----------------------------------------------------------------------
!! *** ROUTINE eos_init ***
!!
!! ** Purpose : initializations for the equation of state
!!
!! ** Method : Read the namelist nameos and control the parameters
!!----------------------------------------------------------------------
INTEGER :: ios ! local integer
INTEGER :: ioptio ! local integer
!!
NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, ln_LEOS, rn_a0, rn_b0, &
& rn_lambda1, rn_mu1, rn_lambda2, rn_mu2, rn_nu
!!----------------------------------------------------------------------
!
READ ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in reference namelist' )
!
READ ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nameos in configuration namelist' )
IF(lwm) WRITE( numond, nameos )
!
rho0 = 1027.51_wp !: volumic mass of reference [kg/m3]
rcp = 3974.00_wp !: heat capacity [J/K]
!
IF(lwp) THEN ! Control print
WRITE(numout,*)
WRITE(numout,*) 'eos_init : equation of state'
WRITE(numout,*) '~~~~~~~~'
WRITE(numout,*) ' Namelist nameos : Chosen the Equation Of Seawater (EOS)'
WRITE(numout,*) ' TEOS-10 : rho=F(Conservative Temperature, Absolute Salinity, depth) ln_TEOS10 = ', ln_TEOS10
WRITE(numout,*) ' EOS-80 : rho=F(Potential Temperature, Practical Salinity, depth) ln_EOS80 = ', ln_EOS80
WRITE(numout,*) ' S-EOS : rho=F(Conservative Temperature, Absolute Salinity, depth) ln_SEOS = ', ln_SEOS
WRITE(numout,*) ' L-EOS : rho=F(Potential Temperature, Practical Salinity, depth) ln_LEOS = ', ln_LEOS
ENDIF
! Check options for equation of state & set neos based on logical flags
ioptio = 0
IF( ln_TEOS10 ) THEN ; ioptio = ioptio+1 ; neos = np_teos10 ; ENDIF
IF( ln_EOS80 ) THEN ; ioptio = ioptio+1 ; neos = np_eos80 ; ENDIF
IF( ln_SEOS ) THEN ; ioptio = ioptio+1 ; neos = np_seos ; ENDIF
IF( ln_LEOS ) THEN ; ioptio = ioptio+1 ; neos = np_leos ; ENDIF
IF( ioptio /= 1 ) CALL ctl_stop("Exactly one equation of state option must be selected")
!
SELECT CASE( neos ) ! check option
!
CASE( np_teos10 ) !== polynomial TEOS-10 ==!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>>> use of TEOS-10 equation of state (cons. temp. and abs. salinity)'
!
l_useCT = .TRUE. ! model temperature is Conservative temperature
!
rdeltaS = 32._wp
r1_S0 = 0.875_wp/35.16504_wp
r1_T0 = 1._wp/40._wp
r1_Z0 = 1.e-4_wp
!
EOS000 = 8.0189615746e+02_wp
EOS100 = 8.6672408165e+02_wp
EOS200 = -1.7864682637e+03_wp
EOS300 = 2.0375295546e+03_wp
EOS400 = -1.2849161071e+03_wp
EOS500 = 4.3227585684e+02_wp
EOS600 = -6.0579916612e+01_wp
EOS010 = 2.6010145068e+01_wp
EOS110 = -6.5281885265e+01_wp
EOS210 = 8.1770425108e+01_wp
EOS310 = -5.6888046321e+01_wp
EOS410 = 1.7681814114e+01_wp
EOS510 = -1.9193502195_wp
EOS020 = -3.7074170417e+01_wp
EOS120 = 6.1548258127e+01_wp
EOS220 = -6.0362551501e+01_wp
EOS320 = 2.9130021253e+01_wp
EOS420 = -5.4723692739_wp
EOS030 = 2.1661789529e+01_wp
EOS130 = -3.3449108469e+01_wp
EOS230 = 1.9717078466e+01_wp
EOS330 = -3.1742946532_wp
EOS040 = -8.3627885467_wp
EOS140 = 1.1311538584e+01_wp
EOS240 = -5.3563304045_wp
EOS050 = 5.4048723791e-01_wp
EOS150 = 4.8169980163e-01_wp
EOS060 = -1.9083568888e-01_wp
EOS001 = 1.9681925209e+01_wp
EOS101 = -4.2549998214e+01_wp
EOS201 = 5.0774768218e+01_wp
EOS301 = -3.0938076334e+01_wp
EOS401 = 6.6051753097_wp
EOS011 = -1.3336301113e+01_wp
EOS111 = -4.4870114575_wp
EOS211 = 5.0042598061_wp
EOS311 = -6.5399043664e-01_wp
EOS021 = 6.7080479603_wp
EOS121 = 3.5063081279_wp
EOS221 = -1.8795372996_wp
EOS031 = -2.4649669534_wp
EOS131 = -5.5077101279e-01_wp
EOS041 = 5.5927935970e-01_wp
EOS002 = 2.0660924175_wp
EOS102 = -4.9527603989_wp
EOS202 = 2.5019633244_wp
EOS012 = 2.0564311499_wp
EOS112 = -2.1311365518e-01_wp
EOS022 = -1.2419983026_wp
EOS003 = -2.3342758797e-02_wp
EOS103 = -1.8507636718e-02_wp
EOS013 = 3.7969820455e-01_wp
!
ALP000 = -6.5025362670e-01_wp
ALP100 = 1.6320471316_wp
ALP200 = -2.0442606277_wp
ALP300 = 1.4222011580_wp
ALP400 = -4.4204535284e-01_wp
ALP500 = 4.7983755487e-02_wp
ALP010 = 1.8537085209_wp
ALP110 = -3.0774129064_wp
ALP210 = 3.0181275751_wp
ALP310 = -1.4565010626_wp
ALP410 = 2.7361846370e-01_wp
ALP020 = -1.6246342147_wp
ALP120 = 2.5086831352_wp
ALP220 = -1.4787808849_wp
ALP320 = 2.3807209899e-01_wp
ALP030 = 8.3627885467e-01_wp
ALP130 = -1.1311538584_wp
ALP230 = 5.3563304045e-01_wp
ALP040 = -6.7560904739e-02_wp
ALP140 = -6.0212475204e-02_wp
ALP050 = 2.8625353333e-02_wp
ALP001 = 3.3340752782e-01_wp
ALP101 = 1.1217528644e-01_wp
ALP201 = -1.2510649515e-01_wp
ALP301 = 1.6349760916e-02_wp
ALP011 = -3.3540239802e-01_wp
ALP111 = -1.7531540640e-01_wp
ALP211 = 9.3976864981e-02_wp
ALP021 = 1.8487252150e-01_wp
ALP121 = 4.1307825959e-02_wp
ALP031 = -5.5927935970e-02_wp
ALP002 = -5.1410778748e-02_wp
ALP102 = 5.3278413794e-03_wp
ALP012 = 6.2099915132e-02_wp
ALP003 = -9.4924551138e-03_wp
!
BET000 = 1.0783203594e+01_wp
BET100 = -4.4452095908e+01_wp
BET200 = 7.6048755820e+01_wp
BET300 = -6.3944280668e+01_wp
BET400 = 2.6890441098e+01_wp
BET500 = -4.5221697773_wp
BET010 = -8.1219372432e-01_wp
BET110 = 2.0346663041_wp
BET210 = -2.1232895170_wp
BET310 = 8.7994140485e-01_wp
BET410 = -1.1939638360e-01_wp
BET020 = 7.6574242289e-01_wp
BET120 = -1.5019813020_wp
BET220 = 1.0872489522_wp
BET320 = -2.7233429080e-01_wp
BET030 = -4.1615152308e-01_wp
BET130 = 4.9061350869e-01_wp
BET230 = -1.1847737788e-01_wp
BET040 = 1.4073062708e-01_wp
BET140 = -1.3327978879e-01_wp
BET050 = 5.9929880134e-03_wp
BET001 = -5.2937873009e-01_wp
BET101 = 1.2634116779_wp
BET201 = -1.1547328025_wp
BET301 = 3.2870876279e-01_wp
BET011 = -5.5824407214e-02_wp
BET111 = 1.2451933313e-01_wp
BET211 = -2.4409539932e-02_wp
BET021 = 4.3623149752e-02_wp
BET121 = -4.6767901790e-02_wp
BET031 = -6.8523260060e-03_wp
BET002 = -6.1618945251e-02_wp
BET102 = 6.2255521644e-02_wp
BET012 = -2.6514181169e-03_wp
BET003 = -2.3025968587e-04_wp
!
PEN000 = -9.8409626043_wp
PEN100 = 2.1274999107e+01_wp
PEN200 = -2.5387384109e+01_wp
PEN300 = 1.5469038167e+01_wp
PEN400 = -3.3025876549_wp
PEN010 = 6.6681505563_wp
PEN110 = 2.2435057288_wp
PEN210 = -2.5021299030_wp
PEN310 = 3.2699521832e-01_wp
PEN020 = -3.3540239802_wp
PEN120 = -1.7531540640_wp
PEN220 = 9.3976864981e-01_wp
PEN030 = 1.2324834767_wp
PEN130 = 2.7538550639e-01_wp
PEN040 = -2.7963967985e-01_wp
PEN001 = -1.3773949450_wp
PEN101 = 3.3018402659_wp
PEN201 = -1.6679755496_wp
PEN011 = -1.3709540999_wp
PEN111 = 1.4207577012e-01_wp
PEN021 = 8.2799886843e-01_wp
PEN002 = 1.7507069098e-02_wp
PEN102 = 1.3880727538e-02_wp
PEN012 = -2.8477365341e-01_wp
!
APE000 = -1.6670376391e-01_wp
APE100 = -5.6087643219e-02_wp
APE200 = 6.2553247576e-02_wp
APE300 = -8.1748804580e-03_wp
APE010 = 1.6770119901e-01_wp
APE110 = 8.7657703198e-02_wp
APE210 = -4.6988432490e-02_wp
APE020 = -9.2436260751e-02_wp
APE120 = -2.0653912979e-02_wp
APE030 = 2.7963967985e-02_wp
APE001 = 3.4273852498e-02_wp
APE101 = -3.5518942529e-03_wp
APE011 = -4.1399943421e-02_wp
APE002 = 7.1193413354e-03_wp
!
BPE000 = 2.6468936504e-01_wp
BPE100 = -6.3170583896e-01_wp
BPE200 = 5.7736640125e-01_wp
BPE300 = -1.6435438140e-01_wp
BPE010 = 2.7912203607e-02_wp
BPE110 = -6.2259666565e-02_wp
BPE210 = 1.2204769966e-02_wp
BPE020 = -2.1811574876e-02_wp
BPE120 = 2.3383950895e-02_wp
BPE030 = 3.4261630030e-03_wp
BPE001 = 4.1079296834e-02_wp
BPE101 = -4.1503681096e-02_wp
BPE011 = 1.7676120780e-03_wp
BPE002 = 1.7269476440e-04_wp
!
CASE( np_eos80 ) !== polynomial EOS-80 formulation ==!
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>>> use of EOS-80 equation of state (pot. temp. and pract. salinity)'
!
l_useCT = .FALSE. ! model temperature is Potential temperature
rdeltaS = 20._wp
r1_S0 = 1._wp/40._wp
r1_T0 = 1._wp/40._wp
r1_Z0 = 1.e-4_wp
!
EOS000 = 9.5356891948e+02_wp
EOS100 = 1.7136499189e+02_wp
EOS200 = -3.7501039454e+02_wp
EOS300 = 5.1856810420e+02_wp
EOS400 = -3.7264470465e+02_wp
EOS500 = 1.4302533998e+02_wp
EOS600 = -2.2856621162e+01_wp
EOS010 = 1.0087518651e+01_wp
EOS110 = -1.3647741861e+01_wp
EOS210 = 8.8478359933_wp
EOS310 = -7.2329388377_wp
EOS410 = 1.4774410611_wp
EOS510 = 2.0036720553e-01_wp
EOS020 = -2.5579830599e+01_wp
EOS120 = 2.4043512327e+01_wp
EOS220 = -1.6807503990e+01_wp
EOS320 = 8.3811577084_wp
EOS420 = -1.9771060192_wp
EOS030 = 1.6846451198e+01_wp
EOS130 = -2.1482926901e+01_wp
EOS230 = 1.0108954054e+01_wp
EOS330 = -6.2675951440e-01_wp
EOS040 = -8.0812310102_wp
EOS140 = 1.0102374985e+01_wp
EOS240 = -4.8340368631_wp
EOS050 = 1.2079167803_wp
EOS150 = 1.1515380987e-01_wp
EOS060 = -2.4520288837e-01_wp
EOS001 = 1.0748601068e+01_wp
EOS101 = -1.7817043500e+01_wp
EOS201 = 2.2181366768e+01_wp
EOS301 = -1.6750916338e+01_wp
EOS401 = 4.1202230403_wp
EOS011 = -1.5852644587e+01_wp
EOS111 = -7.6639383522e-01_wp
EOS211 = 4.1144627302_wp
EOS311 = -6.6955877448e-01_wp
EOS021 = 9.9994861860_wp
EOS121 = -1.9467067787e-01_wp
EOS221 = -1.2177554330_wp
EOS031 = -3.4866102017_wp
EOS131 = 2.2229155620e-01_wp
EOS041 = 5.9503008642e-01_wp
EOS002 = 1.0375676547_wp
EOS102 = -3.4249470629_wp
EOS202 = 2.0542026429_wp
EOS012 = 2.1836324814_wp
EOS112 = -3.4453674320e-01_wp
EOS022 = -1.2548163097_wp
EOS003 = 1.8729078427e-02_wp
EOS103 = -5.7238495240e-02_wp
EOS013 = 3.8306136687e-01_wp
!
ALP000 = -2.5218796628e-01_wp
ALP100 = 3.4119354654e-01_wp
ALP200 = -2.2119589983e-01_wp
ALP300 = 1.8082347094e-01_wp
ALP400 = -3.6936026529e-02_wp
ALP500 = -5.0091801383e-03_wp
ALP010 = 1.2789915300_wp
ALP110 = -1.2021756164_wp
ALP210 = 8.4037519952e-01_wp
ALP310 = -4.1905788542e-01_wp
ALP410 = 9.8855300959e-02_wp
ALP020 = -1.2634838399_wp
ALP120 = 1.6112195176_wp
ALP220 = -7.5817155402e-01_wp
ALP320 = 4.7006963580e-02_wp
ALP030 = 8.0812310102e-01_wp
ALP130 = -1.0102374985_wp
ALP230 = 4.8340368631e-01_wp
ALP040 = -1.5098959754e-01_wp
ALP140 = -1.4394226233e-02_wp
ALP050 = 3.6780433255e-02_wp
ALP001 = 3.9631611467e-01_wp
ALP101 = 1.9159845880e-02_wp
ALP201 = -1.0286156825e-01_wp
ALP301 = 1.6738969362e-02_wp
ALP011 = -4.9997430930e-01_wp
ALP111 = 9.7335338937e-03_wp
ALP211 = 6.0887771651e-02_wp
ALP021 = 2.6149576513e-01_wp
ALP121 = -1.6671866715e-02_wp
ALP031 = -5.9503008642e-02_wp
ALP002 = -5.4590812035e-02_wp
ALP102 = 8.6134185799e-03_wp
ALP012 = 6.2740815484e-02_wp
ALP003 = -9.5765341718e-03_wp
!
BET000 = 2.1420623987_wp
BET100 = -9.3752598635_wp
BET200 = 1.9446303907e+01_wp
BET300 = -1.8632235232e+01_wp
BET400 = 8.9390837485_wp
BET500 = -1.7142465871_wp
BET010 = -1.7059677327e-01_wp
BET110 = 2.2119589983e-01_wp
BET210 = -2.7123520642e-01_wp
BET310 = 7.3872053057e-02_wp
BET410 = 1.2522950346e-02_wp
BET020 = 3.0054390409e-01_wp
BET120 = -4.2018759976e-01_wp
BET220 = 3.1429341406e-01_wp
BET320 = -9.8855300959e-02_wp
BET030 = -2.6853658626e-01_wp
BET130 = 2.5272385134e-01_wp
BET230 = -2.3503481790e-02_wp
BET040 = 1.2627968731e-01_wp
BET140 = -1.2085092158e-01_wp
BET050 = 1.4394226233e-03_wp
BET001 = -2.2271304375e-01_wp
BET101 = 5.5453416919e-01_wp
BET201 = -6.2815936268e-01_wp
BET301 = 2.0601115202e-01_wp
BET011 = -9.5799229402e-03_wp
BET111 = 1.0286156825e-01_wp
BET211 = -2.5108454043e-02_wp
BET021 = -2.4333834734e-03_wp
BET121 = -3.0443885826e-02_wp
BET031 = 2.7786444526e-03_wp
BET002 = -4.2811838287e-02_wp
BET102 = 5.1355066072e-02_wp
BET012 = -4.3067092900e-03_wp
BET003 = -7.1548119050e-04_wp
!
PEN000 = -5.3743005340_wp
PEN100 = 8.9085217499_wp
PEN200 = -1.1090683384e+01_wp
PEN300 = 8.3754581690_wp
PEN400 = -2.0601115202_wp
PEN010 = 7.9263222935_wp
PEN110 = 3.8319691761e-01_wp
PEN210 = -2.0572313651_wp
PEN310 = 3.3477938724e-01_wp
PEN020 = -4.9997430930_wp
PEN120 = 9.7335338937e-02_wp
PEN220 = 6.0887771651e-01_wp
PEN030 = 1.7433051009_wp
PEN130 = -1.1114577810e-01_wp
PEN040 = -2.9751504321e-01_wp
PEN001 = -6.9171176978e-01_wp
PEN101 = 2.2832980419_wp
PEN201 = -1.3694684286_wp
PEN011 = -1.4557549876_wp
PEN111 = 2.2969116213e-01_wp
PEN021 = 8.3654420645e-01_wp
PEN002 = -1.4046808820e-02_wp
PEN102 = 4.2928871430e-02_wp
PEN012 = -2.8729602515e-01_wp
!
APE000 = -1.9815805734e-01_wp
APE100 = -9.5799229402e-03_wp
APE200 = 5.1430784127e-02_wp
APE300 = -8.3694846809e-03_wp
APE010 = 2.4998715465e-01_wp
APE110 = -4.8667669469e-03_wp
APE210 = -3.0443885826e-02_wp
APE020 = -1.3074788257e-01_wp
APE120 = 8.3359333577e-03_wp
APE030 = 2.9751504321e-02_wp
APE001 = 3.6393874690e-02_wp
APE101 = -5.7422790533e-03_wp
APE011 = -4.1827210323e-02_wp
APE002 = 7.1824006288e-03_wp
!
BPE000 = 1.1135652187e-01_wp
BPE100 = -2.7726708459e-01_wp
BPE200 = 3.1407968134e-01_wp
BPE300 = -1.0300557601e-01_wp
BPE010 = 4.7899614701e-03_wp
BPE110 = -5.1430784127e-02_wp
BPE210 = 1.2554227021e-02_wp
BPE020 = 1.2166917367e-03_wp
BPE120 = 1.5221942913e-02_wp
BPE030 = -1.3893222263e-03_wp
BPE001 = 2.8541225524e-02_wp
BPE101 = -3.4236710714e-02_wp
BPE011 = 2.8711395266e-03_wp
BPE002 = 5.3661089288e-04_wp
!
CASE( np_seos ) !== Simplified EOS ==!
r1_S0 = 0.875_wp/35.16504_wp ! Used to convert CT in potential temperature when using bulk formulae (eos_pt_from_ct)
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) ' ==>>> use of simplified eos: '
WRITE(numout,*) ' rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT '
WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rho0'
WRITE(numout,*) ' with the following coefficients :'
WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0
WRITE(numout,*) ' saline cont. coef. rn_b0 = ', rn_b0
WRITE(numout,*) ' cabbeling coef. rn_lambda1 = ', rn_lambda1
WRITE(numout,*) ' cabbeling coef. rn_lambda2 = ', rn_lambda2
WRITE(numout,*) ' thermobar. coef. rn_mu1 = ', rn_mu1
WRITE(numout,*) ' thermobar. coef. rn_mu2 = ', rn_mu2
WRITE(numout,*) ' 2nd cabbel. coef. rn_nu = ', rn_nu
WRITE(numout,*) ' Caution: rn_beta0=0 incompatible with ddm parameterization '
ENDIF
l_useCT = .TRUE. ! Use conservative temperature
!
CASE( np_leos ) !== Linear ISOMIP EOS ==!
r1_S0 = 0.875_wp/35.16504_wp ! Used to convert CT in potential temperature when using bulk formulae (eos_pt_from_ct)
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) ' use of linear ISOMIP eos: rhd(dT=T-(-1),dS=S-(34.2),Z) = '
WRITE(numout,*) ' [ -a0*dT + b0*dS ]/rho0'
WRITE(numout,*)
WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0
WRITE(numout,*) ' saline cont. coef. rn_b0 = ', rn_b0
ENDIF
l_useCT = .TRUE. ! Use conservative temperature
!
CASE DEFAULT !== ERROR in neos ==!
WRITE(ctmp1,*) ' bad flag value for neos = ', neos, '. You should never see this error'
CALL ctl_stop( ctmp1 )
!
END SELECT
!
rho0_rcp = rho0 * rcp
r1_rho0 = 1._wp / rho0
r1_rcp = 1._wp / rcp
r1_rho0_rcp = 1._wp / rho0_rcp
!
IF(lwp) THEN
IF( l_useCT ) THEN
WRITE(numout,*)
WRITE(numout,*) ' ==>>> model uses Conservative Temperature'
WRITE(numout,*) ' Important: model must be initialized with CT and SA fields'
ELSE
WRITE(numout,*)
WRITE(numout,*) ' ==>>> model does not use Conservative Temperature'
ENDIF
ENDIF
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' Associated physical constant'
IF(lwp) WRITE(numout,*) ' volumic mass of reference rho0 = ', rho0 , ' kg/m^3'
IF(lwp) WRITE(numout,*) ' 1. / rho0 r1_rho0 = ', r1_rho0, ' m^3/kg'
IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin'
IF(lwp) WRITE(numout,*) ' rho0 * rcp rho0_rcp = ', rho0_rcp
IF(lwp) WRITE(numout,*) ' 1. / ( rho0 * rcp ) r1_rho0_rcp = ', r1_rho0_rcp
!
END SUBROUTINE eos_init
!!======================================================================
END MODULE eosbn2
tests/ISOMIP+/MY_SRC/isf_oce.F90
View file @ 4dbfc9bb
......@@ -65,7 +65,7 @@ MODULE isf_oce
!
REAL(wp), PARAMETER, PUBLIC :: rLfusisf = 0.334e6_wp !: latent heat of fusion of ice shelf [J/kg]
REAL(wp), PARAMETER, PUBLIC :: rcpisf = 2000.0_wp !: specific heat of ice shelf [J/kg/K]
REAL(wp), PARAMETER, PUBLIC :: rkappa = 0.0_wp !: ISOMIP+ no heat diffusivity through the ice-shelf [m2/s]
REAL(wp), PARAMETER, PUBLIC :: rkappa = 0._wp !: ISOMIP: no heat diffusivity [m2/s]
REAL(wp), PARAMETER, PUBLIC :: rhoisf = 920.0_wp !: volumic mass of ice shelf [kg/m3]
REAL(wp), PARAMETER, PUBLIC :: rtsurf = -20.0 !: surface temperature [C]
!
......@@ -80,8 +80,8 @@ MODULE isf_oce
!
! 2.2 -------- ice shelf cavity melt namelist parameter -------------
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mskisf_cav !:
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_cav , misfkb_cav !:
REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_cav, rfrac_tbl_cav !:
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_cav , misfkb_cav !: shallowest and deepest level of the ice shelf
REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_cav, rfrac_tbl_cav !: thickness and fraction of deepest cell affected by the ice shelf
REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_cav , fwfisf_cav_b !: before and now net fwf from the ice shelf [kg/m2/s]
REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_cav_tsc , risf_cav_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]
TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfcav_fwf !:
......@@ -237,17 +237,28 @@ CONTAINS
!
ierr = 0 ! set to zero if no array to be allocated
!
ALLOCATE( fwfisf_par (jpi,jpj) , fwfisf_par_b(jpi,jpj) , &
& fwfisf_cav (jpi,jpj) , fwfisf_cav_b(jpi,jpj) , &
ALLOCATE( fwfisf_par (jpi,jpj) , fwfisf_cav (jpi,jpj) , &
& fwfisf_oasis(jpi,jpj) , STAT=ialloc )
ierr = ierr + ialloc
!
ALLOCATE( risf_par_tsc(jpi,jpj,jpts) , risf_par_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
ALLOCATE( risf_par_tsc(jpi,jpj,jpts) , STAT=ialloc )
ierr = ierr + ialloc
!
ALLOCATE( risf_cav_tsc(jpi,jpj,jpts) , risf_cav_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
ALLOCATE( risf_cav_tsc(jpi,jpj,jpts) , STAT=ialloc )
ierr = ierr + ialloc
!
#if ! defined key_RK3
! MLF : need to allocate before arrays
ALLOCATE( fwfisf_par_b(jpi,jpj) , fwfisf_cav_b(jpi,jpj) , STAT=ialloc )
ierr = ierr + ialloc
!
ALLOCATE( risf_par_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
ierr = ierr + ialloc
!
ALLOCATE( risf_cav_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
ierr = ierr + ialloc
#endif
!
ALLOCATE( risfload(jpi,jpj) , STAT=ialloc )
ierr = ierr + ialloc
!
......@@ -260,10 +271,17 @@ CONTAINS
! initalisation of fwf and tsc array to 0
risfload (:,:) = 0._wp
fwfisf_oasis(:,:) = 0._wp
#if defined key_RK3
fwfisf_par (:,:) = 0._wp
fwfisf_cav (:,:) = 0._wp
risf_cav_tsc(:,:,:) = 0._wp
risf_par_tsc(:,:,:) = 0._wp
#else
fwfisf_par (:,:) = 0._wp ; fwfisf_par_b (:,:) = 0._wp
fwfisf_cav (:,:) = 0._wp ; fwfisf_cav_b (:,:) = 0._wp
risf_cav_tsc(:,:,:) = 0._wp ; risf_cav_tsc_b(:,:,:) = 0._wp
risf_par_tsc(:,:,:) = 0._wp ; risf_par_tsc_b(:,:,:) = 0._wp
#endif
!
END SUBROUTINE isf_alloc
......
tests/ISOMIP+/MY_SRC/isfcavgam.F90 deleted 100644 → 0
View file @ 0b9e2937
MODULE isfcavgam
!!======================================================================
!! *** MODULE isfgammats ***
!! Ice shelf gamma module : compute exchange coeficient at the ice/ocean interface
!!======================================================================
!! History : 4.1 ! (P. Mathiot) original
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! isfcav_gammats : compute exchange coeficient gamma
!!----------------------------------------------------------------------
USE isf_oce
USE isfutils, ONLY: debug
USE isftbl , ONLY: isf_tbl
USE oce , ONLY: uu, vv ! ocean dynamics
USE phycst , ONLY: grav, vkarmn ! physical constant
USE eosbn2 , ONLY: eos_rab ! equation of state
USE zdfdrg , ONLY: rCd0_top, r_ke0_top ! vertical physics: top/bottom drag coef.
USE iom , ONLY: iom_put !
USE lib_mpp , ONLY: ctl_stop
USE dom_oce ! ocean space and time domain
USE in_out_manager ! I/O manager
!
IMPLICIT NONE
!
PRIVATE
!
PUBLIC isfcav_gammats
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
!
!!-----------------------------------------------------------------------------------------------------
!! PUBLIC SUBROUTINES
!!-----------------------------------------------------------------------------------------------------
!
SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, pRc, pgt, pgs )
!!----------------------------------------------------------------------
!! ** Purpose : compute the coefficient echange for heat and fwf flux
!!
!! ** Method : select the gamma formulation
!! 3 method available (cst, vel and vel_stab)
!!---------------------------------------------------------------------
!!-------------------------- OUT -------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pgt , pgs ! gamma t and gamma s
!!-------------------------- IN -------------------------------------
INTEGER :: Kmm ! ocean time level index
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pqoce, pqfwf ! isf heat and fwf
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! top boundary layer tracer
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pRc ! Richardson number
!!---------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: zutbl, zvtbl ! top boundary layer velocity
!!---------------------------------------------------------------------
!
!==========================================
! 1.: compute velocity in the tbl if needed
!==========================================
!
SELECT CASE ( cn_gammablk )
CASE ( 'spe' )
! gamma is constant (specified in namelist)
! nothing to do
CASE ('vel', 'vel_stab')
! compute velocity in tbl
CALL isf_tbl(Kmm, uu(:,:,:,Kmm) ,zutbl(:,:),'U', miku, rhisf_tbl_cav)
CALL isf_tbl(Kmm, vv(:,:,:,Kmm) ,zvtbl(:,:),'V', mikv, rhisf_tbl_cav)
!
! mask velocity in tbl with ice shelf mask
zutbl(:,:) = zutbl(:,:) * mskisf_cav(:,:)
zvtbl(:,:) = zvtbl(:,:) * mskisf_cav(:,:)
!
! output
CALL iom_put('utbl',zutbl(:,:))
CALL iom_put('vtbl',zvtbl(:,:))
CASE DEFAULT
CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
END SELECT
!
!==========================================
! 2.: compute gamma
!==========================================
!
SELECT CASE ( cn_gammablk )
CASE ( 'spe' ) ! gamma is constant (specified in namelist)
pgt(:,:) = rn_gammat0
pgs(:,:) = rn_gammas0
CASE ( 'vel' ) ! gamma is proportional to u*
CALL gammats_vel ( zutbl, zvtbl, rCd0_top, r_ke0_top, pgt, pgs )
CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u*
CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pRc, pgt, pgs )
CASE DEFAULT
CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
END SELECT
!
!==========================================
! 3.: output and debug
!==========================================
!
CALL iom_put('isfgammat', pgt(:,:))
CALL iom_put('isfgammas', pgs(:,:))
!
IF (ln_isfdebug) THEN
CALL debug( 'isfcav_gam pgt:', pgt(:,:) )
CALL debug( 'isfcav_gam pgs:', pgs(:,:) )
END IF
!
END SUBROUTINE isfcav_gammats
!
!!-----------------------------------------------------------------------------------------------------
!! PRIVATE SUBROUTINES
!!-----------------------------------------------------------------------------------------------------
!
SUBROUTINE gammats_vel( putbl, pvtbl, pCd, pke2, & ! <<== in
& pgt, pgs ) ! ==>> out gammats [m/s]
!!----------------------------------------------------------------------
!! ** Purpose : compute the coefficient echange coefficient
!!
!! ** Method : gamma is velocity dependent ( gt= gt0 * Ustar )
!!
!! ** Reference : Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016
!!---------------------------------------------------------------------
!!-------------------------- OUT -------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pgt, pgs ! gammat and gammas [m/s]
!!-------------------------- IN -------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: putbl, pvtbl ! velocity in the losch top boundary layer
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pCd ! drag coefficient
REAL(wp), INTENT(in ) :: pke2 ! background velocity
!!---------------------------------------------------------------------
INTEGER :: ji, jj ! loop index
REAL(wp), DIMENSION(jpi,jpj) :: zustar
!!---------------------------------------------------------------------
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! compute ustar (AD15 eq. 27)
zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) ) * mskisf_cav(ji,jj)
!
! Compute gammats
pgt(ji,jj) = zustar(ji,jj) * rn_gammat0
pgs(ji,jj) = zustar(ji,jj) * rn_gammas0
END_2D
!
! output ustar
CALL iom_put('isfustar',zustar(:,:))
!
END SUBROUTINE gammats_vel
SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, pRc, & ! <<== in
& pgt , pgs ) ! ==>> out gammats [m/s]
!!----------------------------------------------------------------------
!! ** Purpose : compute the coefficient echange coefficient
!!
!! ** Method : gamma is velocity dependent and stability dependent
!!
!! ** Reference : Holland and Jenkins, 1999, JPO, p1787-1800
!!---------------------------------------------------------------------
!!-------------------------- OUT -------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pgt, pgs ! gammat and gammas
!!-------------------------- IN -------------------------------------
INTEGER :: Kmm ! ocean time level index
REAL(wp), INTENT(in ) :: pke2 ! background velocity squared
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pqoce, pqfwf ! surface heat flux and fwf flux
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pCd ! drag coeficient
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: putbl, pvtbl ! velocity in the losch top boundary layer
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! tracer in the losch top boundary layer
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pRc ! Richardson number
!!---------------------------------------------------------------------
INTEGER :: ji, jj ! loop index
INTEGER :: ikt ! local integer
REAL(wp) :: zdku, zdkv ! U, V shear
REAL(wp) :: zPr, zSc ! Prandtl and Scmidth number
REAL(wp) :: zmob, zmols ! Monin Obukov length, coriolis factor at T point
REAL(wp) :: zbuofdep, zhnu ! Bouyancy length scale, sublayer tickness
REAL(wp) :: zhmax ! limitation of mol
REAL(wp) :: zetastar ! stability parameter
REAL(wp) :: zgmolet, zgmoles, zgturb ! contribution of modelecular sublayer and turbulence
REAL(wp) :: zcoef ! temporary coef
REAL(wp) :: zdep
REAL(wp) :: zeps = 1.0e-20_wp
REAL(wp), PARAMETER :: zxsiN = 0.052_wp ! dimensionless constant
REAL(wp), PARAMETER :: znu = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
REAL(wp), DIMENSION(2) :: zts, zab
REAL(wp), DIMENSION(jpi,jpj) :: zustar ! friction velocity
!!---------------------------------------------------------------------
!
! compute Pr and Sc number (eq ??)
zPr = 13.8_wp
zSc = 2432.0_wp
!
! compute gamma mole (eq ??)
zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
!
! compute gamma
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
ikt = mikt(ji,jj)
! compute ustar
zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) )
IF( zustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think
pgt(ji,jj) = rn_gammat0
pgs(ji,jj) = rn_gammas0
ELSE
! compute bouyancy
zts(jp_tem) = pttbl(ji,jj)
zts(jp_sal) = pstbl(ji,jj)
zdep = gdepw(ji,jj,ikt,Kmm)
!
CALL eos_rab( zts, zdep, zab, Kmm )
!
! compute length scale (Eq ??)
zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pqfwf(ji,jj) )
!
! compute Monin Obukov Length
! Maximum boundary layer depth (Eq ??)
zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp
!
! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??)
zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
!
! compute eta* (stability parameter) (Eq ??)
zetastar = 1._wp / ( SQRT(1._wp + MAX( 0._wp, zxsiN * zustar(ji,jj) &
& / MAX( 1.e-20, ABS(ff_t(ji,jj)) * zmols * pRc(ji,jj) ) )))
!
! compute the sublayer thickness (Eq ??)
zhnu = 5 * znu / MAX( 1.e-20, zustar(ji,jj) )
!
! compute gamma turb (Eq ??)
zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / MAX( 1.e-10, ABS(ff_t(ji,jj)) * zhnu )) &
& + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
!
! compute gammats
pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
END IF
END_2D
! output ustar
CALL iom_put('isfustar',zustar(:,:))
END SUBROUTINE gammats_vel_stab
END MODULE isfcavgam
tests/ISOMIP+/MY_SRC/isfstp.F90 deleted 100644 → 0
View file @ 0b9e2937
MODULE isfstp
!!======================================================================
!! *** MODULE isfstp ***
!! Ice Shelves : compute iceshelf load, melt and heat flux
!!======================================================================
!! History : 3.2 ! 2011-02 (C.Harris ) Original code isf cav
!! X.X ! 2006-02 (C. Wang ) Original code bg03
!! 3.4 ! 2013-03 (P. Mathiot) Merging + parametrization
!! 4.1 ! 2019-09 (P. Mathiot) Split param/explicit ice shelf and re-organisation
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! isfstp : compute iceshelf melt and heat flux
!!----------------------------------------------------------------------
USE isf_oce ! isf variables
USE isfload, ONLY: isf_load ! ice shelf load
USE isftbl , ONLY: isf_tbl_lvl ! ice shelf boundary layer
USE isfpar , ONLY: isf_par, isf_par_init ! ice shelf parametrisation
USE isfcav , ONLY: isf_cav, isf_cav_init ! ice shelf cavity
USE isfcpl , ONLY: isfcpl_rst_write, isfcpl_init ! isf variables
USE dom_oce ! ocean space and time domain
USE oce , ONLY: ssh ! sea surface height
USE domvvl, ONLY: ln_vvl_zstar ! zstar logical
USE zdfdrg, ONLY: r_Cdmin_top, r_ke0_top ! vertical physics: top/bottom drag coef.
!
USE lib_mpp, ONLY: ctl_stop, ctl_nam
USE fldread, ONLY: FLD, FLD_N
USE in_out_manager ! I/O manager
USE timing
IMPLICIT NONE
PRIVATE
PUBLIC isf_stp, isf_init, isf_nam ! routine called in sbcmod and divhor
!! * Substitutions
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: isfstp.F90 15574 2021-12-03 19:32:50Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE isf_stp( kt, Kmm )
!!---------------------------------------------------------------------
!! *** ROUTINE isf_stp ***
!!
!! ** Purpose : compute total heat flux and total fwf due to ice shelf melt
!!
!! ** Method : For each case (parametrisation or explicity cavity) :
!! - define the before fields
!! - compute top boundary layer properties
!! (in case of parametrisation, this is the
!! depth range model array used to compute mean far fields properties)
!! - compute fluxes
!! - write restart variables
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time step
INTEGER, INTENT(in) :: Kmm ! ocean time level index
!
INTEGER :: jk ! loop index
#if defined key_qco
REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t ! 3D workspace
#endif
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('isf')
!
!=======================================================================
! 1.: compute melt and associated heat fluxes in the ice shelf cavities
!=======================================================================
!
IF ( ln_isfcav_mlt ) THEN
!
! 1.1: before time step
IF ( kt /= nit000 ) THEN
risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:)
fwfisf_cav_b(:,:) = fwfisf_cav(:,:)
END IF
!
! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:)
#if defined key_qco
DO jk = 1, jpk
ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
END DO
CALL isf_tbl_lvl( ht(:,:), ze3t , misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav )
#else
CALL isf_tbl_lvl( ht(:,:), e3t(:,:,:,Kmm), misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav )
#endif
!
! 1.3: compute ice shelf melt
CALL isf_cav( kt, Kmm, risf_cav_tsc, fwfisf_cav )
!
END IF
!
!=================================================================================
! 2.: compute melt and associated heat fluxes for not resolved ice shelf cavities
!=================================================================================
!
IF ( ln_isfpar_mlt ) THEN
!
! 2.1: before time step
IF ( kt /= nit000 ) THEN
risf_par_tsc_b(:,:,:) = risf_par_tsc(:,:,:)
fwfisf_par_b (:,:) = fwfisf_par (:,:)
END IF
!
! 2.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
! by simplicity, we assume the top level where param applied do not change with time (done in init part)
rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
#if defined key_qco
DO jk = 1, jpk
ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
END DO
CALL isf_tbl_lvl( ht(:,:), ze3t , misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
#else
CALL isf_tbl_lvl( ht(:,:), e3t(:,:,:,Kmm), misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
#endif
!
! 2.3: compute ice shelf melt
CALL isf_par( kt, Kmm, risf_par_tsc, fwfisf_par )
!
END IF
!
!==================================================================================
! 3.: output specific restart variable in case of coupling with an ice sheet model
!==================================================================================
!
IF ( ln_isfcpl .AND. lrst_oce ) CALL isfcpl_rst_write(kt, Kmm)
!
IF( ln_timing ) CALL timing_stop('isf')
!
END SUBROUTINE isf_stp
SUBROUTINE isf_init( Kbb, Kmm, Kaa )
!!---------------------------------------------------------------------
!! *** ROUTINE isfstp_init ***
!!
!! ** Purpose : Initialisation of the ice shelf public variables
!!
!! ** Method : Read the namisf namelist, check option compatibility and set derived parameters
!!
!! ** Action : - read namisf parameters
!! - allocate memory
!! - output print
!! - ckeck option compatibility
!! - call cav/param/isfcpl init routine
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
!!----------------------------------------------------------------------
!
! constrain: l_isfoasis need to be known
!
CALL isf_nam() ! Read namelist
!
CALL isf_alloc() ! Allocate public array
!
CALL isf_ctl() ! check option compatibility
!
IF( ln_isfcav ) CALL isf_load( Kmm, risfload ) ! compute ice shelf load
!
! terminate routine now if no ice shelf melt formulation specify
IF( ln_isf ) THEN
!
IF( ln_isfcav_mlt ) CALL isf_cav_init() ! initialisation melt in the cavity
!
IF( ln_isfpar_mlt ) CALL isf_par_init() ! initialisation parametrised melt
!
IF( ln_isfcpl ) CALL isfcpl_init( Kbb, Kmm, Kaa ) ! initialisation ice sheet coupling
!
END IF
END SUBROUTINE isf_init
SUBROUTINE isf_ctl()
!!---------------------------------------------------------------------
!! *** ROUTINE isf_ctl ***
!!
!! ** Purpose : output print and option compatibility check
!!
!!----------------------------------------------------------------------
IF (lwp) THEN
WRITE(numout,*)
WRITE(numout,*) 'isf_init : ice shelf initialisation'
WRITE(numout,*) '~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namisf :'
!
WRITE(numout,*) ' ice shelf cavity (open or parametrised) ln_isf = ', ln_isf
WRITE(numout,*)
!
IF ( ln_isf ) THEN
#if key_qco
# if ! defined key_isf
CALL ctl_stop( 'STOP', 'isf_ctl: ice shelf requires both ln_isf=T AND key_isf activated' )
# endif
#endif
WRITE(numout,*) ' Add debug print in isf module ln_isfdebug = ', ln_isfdebug
WRITE(numout,*)
WRITE(numout,*) ' melt inside the cavity ln_isfcav_mlt = ', ln_isfcav_mlt
IF ( ln_isfcav_mlt) THEN
WRITE(numout,*) ' melt formulation cn_isfcav_mlt= ', TRIM(cn_isfcav_mlt)
WRITE(numout,*) ' thickness of the top boundary layer rn_htbl = ', rn_htbl
WRITE(numout,*) ' gamma formulation cn_gammablk = ', TRIM(cn_gammablk)
IF ( TRIM(cn_gammablk) .NE. 'spe' ) THEN
WRITE(numout,*) ' gammat coefficient rn_gammat0 = ', rn_gammat0
WRITE(numout,*) ' gammas coefficient rn_gammas0 = ', rn_gammas0
WRITE(numout,*) ' top background ke used (from namdrg_top) rn_ke0 = ', r_ke0_top
WRITE(numout,*) ' top drag coef. used (from namdrg_top) rn_Cd0 = ', r_Cdmin_top
END IF
END IF
WRITE(numout,*) ''
!
WRITE(numout,*) ' ice shelf melt parametrisation ln_isfpar_mlt = ', ln_isfpar_mlt
IF ( ln_isfpar_mlt ) THEN
WRITE(numout,*) ' isf parametrisation formulation cn_isfpar_mlt = ', TRIM(cn_isfpar_mlt)
END IF
WRITE(numout,*) ''
!
WRITE(numout,*) ' Coupling to an ice sheet model ln_isfcpl = ', ln_isfcpl
IF ( ln_isfcpl ) THEN
WRITE(numout,*) ' conservation activated ln_isfcpl_cons = ', ln_isfcpl_cons
WRITE(numout,*) ' number of call of the extrapolation loop = ', nn_drown
ENDIF
WRITE(numout,*) ''
!
ELSE
!
IF ( ln_isfcav ) THEN
WRITE(numout,*) ''
WRITE(numout,*) ' W A R N I N G: ice shelf cavities are open BUT no melt will be computed or read from file !'
WRITE(numout,*) ''
END IF
!
END IF
IF (ln_isfcav) THEN
WRITE(numout,*) ' Ice shelf load method cn_isfload = ', TRIM(cn_isfload)
WRITE(numout,*) ' Temperature used to compute the ice shelf load = ', rn_isfload_T
WRITE(numout,*) ' Salinity used to compute the ice shelf load = ', rn_isfload_S
END IF
WRITE(numout,*) ''
FLUSH(numout)
END IF
!
!---------------------------------------------------------------------------------------------------------------------
! sanity check ! issue ln_isfcav not yet known as well as l_isfoasis => move this call in isf_stp ?
! melt in the cavity without cavity
IF ( ln_isfcav_mlt .AND. (.NOT. ln_isfcav) ) &
& CALL ctl_stop('ice shelf melt in the cavity activated (ln_isfcav_mlt) but no cavity detected in domcfg (ln_isfcav), STOP' )
!
! ice sheet coupling without cavity
IF ( ln_isfcpl .AND. (.NOT. ln_isfcav) ) &
& CALL ctl_stop('coupling with an ice sheet model detected (ln_isfcpl) but no cavity detected in domcfg (ln_isfcav), STOP' )
!
IF ( ln_isfcpl .AND. ln_isfcpl_cons .AND. ln_linssh ) &
& CALL ctl_stop( 'The coupling between NEMO and an ice sheet model with the conservation option does not work with the linssh option' )
!
IF ( l_isfoasis .AND. .NOT. ln_isf ) CALL ctl_stop( ' OASIS send ice shelf fluxes to NEMO but NEMO does not have the isf module activated' )
!
IF ( l_isfoasis .AND. ln_isf ) THEN
!
CALL ctl_stop( 'namelist combination ln_cpl and ln_isf not tested' )
!
! NEMO coupled to ATMO model with isf cavity need oasis method for melt computation
IF ( ln_isfcav_mlt .AND. TRIM(cn_isfcav_mlt) /= 'oasis' ) CALL ctl_stop( 'cn_isfcav_mlt = oasis is the only option availble if fwf send by oasis' )
IF ( ln_isfpar_mlt .AND. TRIM(cn_isfpar_mlt) /= 'oasis' ) CALL ctl_stop( 'cn_isfpar_mlt = oasis is the only option availble if fwf send by oasis' )
!
! oasis melt computation not tested (coded but not tested)
IF ( ln_isfcav_mlt .OR. ln_isfpar_mlt ) THEN
IF ( TRIM(cn_isfcav_mlt) == 'oasis' ) CALL ctl_stop( 'cn_isfcav_mlt = oasis not tested' )
IF ( TRIM(cn_isfpar_mlt) == 'oasis' ) CALL ctl_stop( 'cn_isfpar_mlt = oasis not tested' )
END IF
!
! oasis melt computation with cavity open and cavity parametrised (not coded)
IF ( ln_isfcav_mlt .AND. ln_isfpar_mlt ) THEN
IF ( TRIM(cn_isfpar_mlt) == 'oasis' .AND. TRIM(cn_isfcav_mlt) == 'oasis' ) CALL ctl_stop( 'cn_isfpar_mlt = oasis and cn_isfcav_mlt = oasis not coded' )
END IF
!
! compatibility ice shelf and vvl
IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
!
END IF
END SUBROUTINE isf_ctl
SUBROUTINE isf_nam
!!---------------------------------------------------------------------
!! *** ROUTINE isf_nam ***
!!
!! ** Purpose : Read ice shelf namelist cfg and ref
!!
!!----------------------------------------------------------------------
INTEGER :: ios ! Local integer output status for namelist read
!!----------------------------------------------------------------------
NAMELIST/namisf/ ln_isf , &
& cn_gammablk , rn_gammat0 , rn_gammas0 , rn_htbl, sn_isfcav_fwf, &
& ln_isfcav_mlt , cn_isfcav_mlt , sn_isfcav_fwf , &
& ln_isfpar_mlt , cn_isfpar_mlt , sn_isfpar_fwf , &
& sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff, &
& ln_isfcpl , nn_drown , ln_isfcpl_cons, ln_isfdebug, &
& cn_isfload , rn_isfload_T , rn_isfload_S , cn_isfdir , &
& rn_isfpar_bg03_gt0
!!----------------------------------------------------------------------
!
READ ( numnam_ref, namisf, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namisf in reference namelist' )
!
READ ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namisf in configuration namelist' )
IF(lwm) WRITE ( numond, namisf )
END SUBROUTINE isf_nam
!!
!!======================================================================
END MODULE isfstp
tests/ISOMIP+/MY_SRC/istate.F90
View file @ 4dbfc9bb
......@@ -50,7 +50,7 @@ MODULE istate
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: istate.F90 15581 2021-12-07 13:08:22Z techene $
!! $Id: istate.F90 14991 2021-06-14 19:52:31Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -79,10 +79,12 @@ CONTAINS
CALL dta_tsd_init ! Initialisation of T & S input data
IF( ln_c1d) CALL dta_uvd_init ! Initialisation of U & V input data (c1d only)
rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk
rn2b (:,:,: ) = 0._wp ; rn2 (:,:,: ) = 0._wp ! set one for all to 0 at levels 1 and jpk
ts (:,:,:,:,Kaa) = 0._wp ! set one for all to 0 at level jpk
rab_b(:,:,:,: ) = 0._wp ; rab_n(:,:,:,:) = 0._wp ! set one for all to 0 at level jpk
ts (:,:,:,:,Kaa) = 0._wp ; rn2 (:,:,: ) = 0._wp ! set one for all to 0 at levels 1 and jpk
IF ( ALLOCATED( rhd ) ) THEN ! SWE, for example, will not have allocated these
rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk
rn2b (:,:,: ) = 0._wp ! set one for all to 0 at level jpk
rab_b(:,:,:,: ) = 0._wp ; rab_n(:,:,:,:) = 0._wp ! set one for all to 0 at level jpk
ENDIF
#if defined key_agrif
uu (:,:,: ,Kaa) = 0._wp ! used in agrif_oce_sponge at initialization
vv (:,:,: ,Kaa) = 0._wp ! used in agrif_oce_sponge at initialization
......@@ -113,7 +115,7 @@ CONTAINS
! ! Initialization of ocean to zero
!
IF( ln_tsd_init ) THEN
CALL dta_tsd( nit000, 'ini', ts(:,:,:,:,Kbb) ) ! read 3D T and S data at nit000
CALL dta_tsd( nit000, ts(:,:,:,:,Kbb), 'ini' ) ! read 3D T and S data at nit000
ENDIF
!
IF( ln_uvd_init .AND. ln_c1d ) THEN
......@@ -141,20 +143,6 @@ CONTAINS
ENDIF
#endif
!
! Initialize "now" barotropic velocities:
! Do it whatever the free surface method, these arrays being used eventually
!
!!gm the use of umask & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
#if ! defined key_RK3
uu_b(:,:,Kmm) = 0._wp ; vv_b(:,:,Kmm) = 0._wp
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
END_3D
uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
#endif
!
#if defined key_RK3
IF( .NOT. ln_rstart ) THEN
#endif
......@@ -171,6 +159,25 @@ CONTAINS
!
#if defined key_RK3
ENDIF
#endif
!
! Initialize "now" barotropic velocities:
! Do it whatever the free surface method, these arrays being used eventually
!
#if defined key_RK3
IF( .NOT. ln_rstart ) THEN
uu_b(:,:,Kmm) = uu_b(:,:,Kbb) ! Kmm value set to Kbb for initialisation in Agrif_Regrid in namo_gcm
vv_b(:,:,Kmm) = vv_b(:,:,Kbb)
ENDIF
#else
!!gm the use of umask & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
uu_b(:,:,Kmm) = 0._wp ; vv_b(:,:,Kmm) = 0._wp
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
END_3D
uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
#endif
!
END SUBROUTINE istate_init
......
tests/ISOMIP+/MY_SRC/sbcfwb.F90 deleted 100644 → 0
View file @ 0b9e2937
MODULE sbcfwb
!!======================================================================
!! *** MODULE sbcfwb ***
!! Ocean fluxes : domain averaged freshwater budget
!!======================================================================
!! History : OPA ! 2001-02 (E. Durand) Original code
!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
!! 3.0 ! 2006-08 (G. Madec) Surface module
!! 3.2 ! 2009-07 (C. Talandier) emp mean s spread over erp area
!! 3.6 ! 2014-11 (P. Mathiot ) add ice shelf melting
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! sbc_fwb : freshwater budget for global ocean configurations (free surface & forced mode)
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
USE sbc_oce ! surface ocean boundary condition
USE isf_oce , ONLY : fwfisf_cav, fwfisf_par, ln_isfcpl, ln_isfcpl_cons, risfcpl_cons_ssh ! ice shelf melting contribution
USE sbc_ice , ONLY : snwice_mass, snwice_mass_b, snwice_fmass
USE phycst ! physical constants
USE sbcrnf ! ocean runoffs
USE sbcssr ! Sea-Surface damping terms
!
USE in_out_manager ! I/O manager
USE iom ! IOM
USE lib_mpp ! distribued memory computing library
USE timing ! Timing
USE lbclnk ! ocean lateral boundary conditions
USE lib_fortran !
IMPLICIT NONE
PRIVATE
PUBLIC sbc_fwb ! routine called by step
REAL(wp) :: rn_fwb0 ! initial freshwater adjustment flux [kg/m2/s] (nn_fwb = 2 only)
REAL(wp) :: a_fwb ! annual domain averaged freshwater budget from the previous year
REAL(wp) :: a_fwb_b ! annual domain averaged freshwater budget from the year before or at initial state
REAL(wp) :: a_fwb_ini ! initial domain averaged freshwater budget
REAL(wp) :: area ! global mean ocean surface (interior domain)
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: sbcfwb.F90 11395 2019-08-02 14:19:00Z mathiot $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE sbc_fwb( kt, kn_fwb, kn_fsbc, Kmm )
!!---------------------------------------------------------------------
!! *** ROUTINE sbc_fwb ***
!!
!! ** Purpose : Control the mean sea surface drift
!!
!! ** Method : several ways depending on kn_fwb
!! =0 no control
!! =1 global mean of emp set to zero at each nn_fsbc time step
!! =2 annual global mean corrected from previous year
!! =3 global mean of emp set to zero at each nn_fsbc time step
!! & spread out over erp area depending its sign
!! Note: if sea ice is embedded it is taken into account when computing the budget
!!----------------------------------------------------------------------
INTEGER, INTENT( in ) :: kt ! ocean time-step index
INTEGER, INTENT( in ) :: kn_fsbc !
INTEGER, INTENT( in ) :: kn_fwb ! ocean time-step index
INTEGER, INTENT( in ) :: Kmm ! ocean time level index
!
INTEGER :: ios, inum, ikty ! local integers
REAL(wp) :: z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp ! local scalars
REAL(wp) :: zsurf_neg, zsurf_pos, zsurf_tospread, zcoef ! - -
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmsk_neg, ztmsk_pos, z_wgt ! 2D workspaces
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmsk_tospread, zerp_cor ! - -
REAL(wp) ,DIMENSION(1) :: z_fwfprv
COMPLEX(dp),DIMENSION(1) :: y_fwfnow
!
NAMELIST/namsbc_fwb/rn_fwb0
!!----------------------------------------------------------------------
!
IF( kt == nit000 ) THEN
READ( numnam_ref, namsbc_fwb, IOSTAT = ios, ERR = 901 )
901 IF( ios /= 0 ) CALL ctl_nam( ios, 'namsbc_fwb in reference namelist' )
READ( numnam_cfg, namsbc_fwb, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam( ios, 'namsbc_fwb in configuration namelist' )
IF(lwm) WRITE( numond, namsbc_fwb )
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) 'sbc_fwb : FreshWater Budget correction'
WRITE(numout,*) '~~~~~~~'
IF( kn_fwb == 1 ) WRITE(numout,*) ' instantaneously set to zero'
IF( kn_fwb == 4 ) WRITE(numout,*) ' instantaneously set to zero with heat and salt flux correction (ISOMIP+)'
IF( kn_fwb == 3 ) WRITE(numout,*) ' fwf set to zero and spread out over erp area'
IF( kn_fwb == 2 ) THEN
WRITE(numout,*) ' adjusted from previous year budget'
WRITE(numout,*)
WRITE(numout,*) ' Namelist namsbc_fwb'
WRITE(numout,*) ' Initial freshwater adjustment flux [kg/m2/s] = ', rn_fwb0
END IF
ENDIF
!
IF( kn_fwb == 3 .AND. nn_sssr /= 2 ) CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' )
IF( kn_fwb == 3 .AND. ln_isfcav ) CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' )
!
area = glob_sum( 'sbcfwb', e1e2t(:,:) * tmask(:,:,1)) ! interior global domain surface
! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes
! and in case of no melt, it can generate HSSW.
!
#if ! defined key_si3 && ! defined key_cice
snwice_mass_b(:,:) = 0.e0 ! no sea-ice model is being used : no snow+ice mass
snwice_mass (:,:) = 0.e0
snwice_fmass (:,:) = 0.e0
#endif
!
ENDIF
SELECT CASE ( kn_fwb )
!
CASE ( 1 ) !== global mean fwf set to zero ==!
!
IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
y_fwfnow(1) = local_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) - snwice_fmass(:,:) ) )
CALL mpp_delay_sum( 'sbcfwb', 'fwb', y_fwfnow(:), z_fwfprv(:), kt == nitend - nn_fsbc + 1 )
z_fwfprv(1) = z_fwfprv(1) / area
zcoef = z_fwfprv(1) * rcp
emp(:,:) = emp(:,:) - z_fwfprv(1) * tmask(:,:,1)
qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
! outputs
IF( iom_use('hflx_fwb_cea') ) CALL iom_put( 'hflx_fwb_cea', zcoef * sst_m(:,:) * tmask(:,:,1) )
IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', z_fwfprv(1) * tmask(:,:,1) )
ENDIF
!
CASE ( 2 ) !== fw adjustment based on fw budget at the end of the previous year ==!
! simulation is supposed to start 1st of January
IF( kt == nit000 ) THEN ! initialisation
! ! set the fw adjustment (a_fwb)
IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb_b', ldstop = .FALSE. ) > 0 & ! as read from restart file
& .AND. iom_varid( numror, 'a_fwb', ldstop = .FALSE. ) > 0 ) THEN
IF(lwp) WRITE(numout,*) 'sbc_fwb : reading freshwater-budget from restart file'
CALL iom_get( numror, 'a_fwb_b', a_fwb_b )
CALL iom_get( numror, 'a_fwb' , a_fwb )
!
a_fwb_ini = a_fwb_b
ELSE ! as specified in namelist
IF(lwp) WRITE(numout,*) 'sbc_fwb : setting freshwater-budget from namelist rn_fwb0'
a_fwb = rn_fwb0
a_fwb_b = 0._wp ! used only the first year then it is replaced by a_fwb_ini
!
a_fwb_ini = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) ) &
& * rho0 / ( area * rday * REAL(nyear_len(1), wp) )
END IF
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*)'sbc_fwb : freshwater-budget at the end of previous year = ', a_fwb , 'kg/m2/s'
IF(lwp) WRITE(numout,*)' freshwater-budget at initial state = ', a_fwb_ini, 'kg/m2/s'
!
ELSE
! at the end of year n:
ikty = nyear_len(1) * 86400 / NINT(rn_Dt)
IF( MOD( kt, ikty ) == 0 ) THEN ! Update a_fwb at the last time step of a year
! It should be the first time step of a year MOD(kt-1,ikty) but then the restart would be wrong
! Hence, we make a small error here but the code is restartable
a_fwb_b = a_fwb_ini
! mean sea level taking into account ice+snow
a_fwb = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) )
a_fwb = a_fwb * rho0 / ( area * rday * REAL(nyear_len(1), wp) ) ! convert in kg/m2/s
ENDIF
!
ENDIF
!
IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! correct the freshwater fluxes using previous year budget minus initial state
zcoef = ( a_fwb - a_fwb_b )
emp(:,:) = emp(:,:) + zcoef * tmask(:,:,1)
qns(:,:) = qns(:,:) - zcoef * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
! outputs
IF( iom_use('hflx_fwb_cea') ) CALL iom_put( 'hflx_fwb_cea', -zcoef * rcp * sst_m(:,:) * tmask(:,:,1) )
IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', -zcoef * tmask(:,:,1) )
ENDIF
! Output restart information
IF( lrst_oce ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'sbc_fwb : writing FW-budget adjustment to ocean restart file at it = ', kt
IF(lwp) WRITE(numout,*) '~~~~'
CALL iom_rstput( kt, nitrst, numrow, 'a_fwb_b', a_fwb_b )
CALL iom_rstput( kt, nitrst, numrow, 'a_fwb', a_fwb )
END IF
!
IF( kt == nitend .AND. lwp ) THEN
WRITE(numout,*) 'sbc_fwb : freshwater-budget at the end of simulation (year now) = ', a_fwb , 'kg/m2/s'
WRITE(numout,*) ' freshwater-budget at initial state = ', a_fwb_b, 'kg/m2/s'
ENDIF
!
CASE ( 3 ) !== global fwf set to zero and spread out over erp area ==!
!
ALLOCATE( ztmsk_neg(jpi,jpj) , ztmsk_pos(jpi,jpj) , ztmsk_tospread(jpi,jpj) , z_wgt(jpi,jpj) , zerp_cor(jpi,jpj) )
!
IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
ztmsk_pos(:,:) = tmask_i(:,:) ! Select <0 and >0 area of erp
WHERE( erp < 0._wp ) ztmsk_pos = 0._wp
ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
! ! fwf global mean (excluding ocean to ice/snow exchanges)
z_fwf = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) - snwice_fmass(:,:) ) ) / area
!
IF( z_fwf < 0._wp ) THEN ! spread out over >0 erp area to increase evaporation
zsurf_pos = glob_sum( 'sbcfwb', e1e2t(:,:)*ztmsk_pos(:,:) )
zsurf_tospread = zsurf_pos
ztmsk_tospread(:,:) = ztmsk_pos(:,:)
ELSE ! spread out over <0 erp area to increase precipitation
zsurf_neg = glob_sum( 'sbcfwb', e1e2t(:,:)*ztmsk_neg(:,:) ) ! Area filled by <0 and >0 erp
zsurf_tospread = zsurf_neg
ztmsk_tospread(:,:) = ztmsk_neg(:,:)
ENDIF
!
zsum_fwf = glob_sum( 'sbcfwb', e1e2t(:,:) * z_fwf ) ! fwf global mean over <0 or >0 erp area
!!gm : zsum_fwf = z_fwf * area ??? it is right? I think so....
z_fwf_nsrf = zsum_fwf / ( zsurf_tospread + rsmall )
! ! weight to respect erp field 2D structure
zsum_erp = glob_sum( 'sbcfwb', ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
! ! final correction term to apply
zerp_cor(:,:) = -1. * z_fwf_nsrf * zsurf_tospread * z_wgt(:,:)
!
!!gm ===>>>> lbc_lnk should be useless as all the computation is done over the whole domain !
CALL lbc_lnk( 'sbcfwb', zerp_cor, 'T', 1.0_wp )
!
emp(:,:) = emp(:,:) + zerp_cor(:,:)
qns(:,:) = qns(:,:) - zerp_cor(:,:) * rcp * sst_m(:,:) ! account for change to the heat budget due to fw correction
erp(:,:) = erp(:,:) + zerp_cor(:,:)
! outputs
IF( iom_use('hflx_fwb_cea') ) CALL iom_put( 'hflx_fwb_cea', -zerp_cor(:,:) * rcp * sst_m(:,:) )
IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', -zerp_cor(:,:) )
!
IF( lwp ) THEN ! control print
IF( z_fwf < 0._wp ) THEN
WRITE(numout,*)' z_fwf < 0'
WRITE(numout,*)' SUM(erp+) = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
ELSE
WRITE(numout,*)' z_fwf >= 0'
WRITE(numout,*)' SUM(erp-) = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
ENDIF
WRITE(numout,*)' SUM(empG) = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv'
WRITE(numout,*)' z_fwf = ', z_fwf ,' Kg/m2/s'
WRITE(numout,*)' z_fwf_nsrf = ', z_fwf_nsrf ,' Kg/m2/s'
WRITE(numout,*)' MIN(zerp_cor) = ', MINVAL(zerp_cor)
WRITE(numout,*)' MAX(zerp_cor) = ', MAXVAL(zerp_cor)
ENDIF
ENDIF
DEALLOCATE( ztmsk_neg , ztmsk_pos , ztmsk_tospread , z_wgt , zerp_cor )
!
CASE ( 4 ) !== global mean fwf set to zero (ISOMIP case) ==!
!
IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
z_fwf = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) - snwice_fmass(:,:) ) )
!
! correction for ice sheet coupling testing (ie remove the excess through the surface)
! test impact on the melt as conservation correction made in depth
! test conservation level as sbcfwb is conserving
! avoid the model to blow up for large ssh drop (isomip OCEAN3 with melt switch off and uniform T/S)
IF (ln_isfcpl .AND. ln_isfcpl_cons) THEN
z_fwf = z_fwf + glob_sum( 'sbcfwb', e1e2t(:,:) * risfcpl_cons_ssh(:,:) * rho0 )
END IF
!
z_fwf = z_fwf / area
zcoef = z_fwf * rcp
emp(:,:) = emp(:,:) - z_fwf * tmask(:,:,1) ! (Eq. 34 AD2015)
qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! (Eq. 35 AD2015) ! use sst_m to avoid generation of any bouyancy fluxes
sfx(:,:) = sfx(:,:) + z_fwf * sss_m(:,:) * tmask(:,:,1) ! (Eq. 36 AD2015) ! use sss_m to avoid generation of any bouyancy fluxes
ENDIF
!
CASE DEFAULT !== you should never be there ==!
CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1, 2 or 3' )
!
END SELECT
!
END SUBROUTINE sbc_fwb
!!======================================================================
END MODULE sbcfwb
tests/ISOMIP+/MY_SRC/tradmp.F90 deleted 100644 → 0
View file @ 0b9e2937
MODULE tradmp
!!======================================================================
!! *** MODULE tradmp ***
!! Ocean physics: internal restoring trend on active tracers (T and S)
!!======================================================================
!! History : OPA ! 1991-03 (O. Marti, G. Madec) Original code
!! ! 1992-06 (M. Imbard) doctor norme
!! ! 1998-07 (M. Imbard, G. Madec) ORCA version
!! 7.0 ! 2001-02 (M. Imbard) add distance to coast, Original code
!! 8.1 ! 2001-02 (G. Madec, E. Durand) cleaning
!! NEMO 1.0 ! 2002-08 (G. Madec, E. Durand) free form + modules
!! 3.2 ! 2009-08 (G. Madec, C. Talandier) DOCTOR norm for namelist parameter
!! 3.3 ! 2010-06 (C. Ethe, G. Madec) merge TRA-TRC
!! 3.4 ! 2011-04 (G. Madec, C. Ethe) Merge of dtatem and dtasal + suppression of CPP keys
!! 3.6 ! 2015-06 (T. Graham) read restoring coefficient in a file
!! 3.7 ! 2015-10 (G. Madec) remove useless trends arrays
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! tra_dmp_alloc : allocate tradmp arrays
!! tra_dmp : update the tracer trend with the internal damping
!! tra_dmp_init : initialization, namlist read, parameters control
!!----------------------------------------------------------------------
USE oce ! ocean: variables
USE dom_oce ! ocean: domain variables
USE trd_oce ! trends: ocean variables
USE trdtra ! trends manager: tracers
USE zdf_oce ! ocean: vertical physics
USE phycst ! physical constants
USE dtatsd ! data: temperature & salinity
USE zdfmxl ! vertical physics: mixed layer depth
!
USE in_out_manager ! I/O manager
USE iom ! XIOS
USE lib_mpp ! MPP library
USE prtctl ! Print control
USE timing ! Timing
IMPLICIT NONE
PRIVATE
PUBLIC tra_dmp ! called by step.F90
PUBLIC tra_dmp_init ! called by nemogcm.F90
! !!* Namelist namtra_dmp : T & S newtonian damping *
LOGICAL , PUBLIC :: ln_tradmp !: internal damping flag
INTEGER , PUBLIC :: nn_zdmp !: = 0/1/2 flag for damping in the mixed layer
CHARACTER(LEN=200) , PUBLIC :: cn_resto !: name of netcdf file containing restoration coefficient field
!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: resto !: restoring coeff. on T and S (s-1)
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: tradmp.F90 15574 2021-12-03 19:32:50Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
INTEGER FUNCTION tra_dmp_alloc()
!!----------------------------------------------------------------------
!! *** FUNCTION tra_dmp_alloc ***
!!----------------------------------------------------------------------
ALLOCATE( resto(jpi,jpj,jpk), STAT= tra_dmp_alloc )
!
CALL mpp_sum ( 'tradmp', tra_dmp_alloc )
IF( tra_dmp_alloc > 0 ) CALL ctl_warn('tra_dmp_alloc: allocation of arrays failed')
!
END FUNCTION tra_dmp_alloc
SUBROUTINE tra_dmp( kt, Kbb, Kmm, pts, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_dmp ***
!!
!! ** Purpose : Compute the tracer trend due to a newtonian damping
!! of the tracer field towards given data field and add it to the
!! general tracer trends.
!!
!! ** Method : Newtonian damping towards t_dta and s_dta computed
!! and add to the general tracer trends:
!! ta = ta + resto * (t_dta - tb)
!! sa = sa + resto * (s_dta - sb)
!! The trend is computed either throughout the water column
!! (nlmdmp=0) or in area of weak vertical mixing (nlmdmp=1) or
!! below the well mixed layer (nlmdmp=2)
!!
!! ** Action : - tsa: tracer trends updated with the damping trend
!!----------------------------------------------------------------------
INTEGER, INTENT(in ) :: kt ! ocean time-step index
INTEGER, INTENT(in ) :: Kbb, Kmm, Krhs ! time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts) :: zts_dta
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zwrk
REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrdts
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('tra_dmp')
!
IF( l_trdtra .OR. iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN !* Save ta and sa trends
ALLOCATE( ztrdts(A2D(nn_hls),jpk,jpts) )
DO jn = 1, jpts
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
ztrdts(ji,jj,jk,jn) = pts(ji,jj,jk,jn,Krhs)
END_3D
END DO
ENDIF
! !== input T-S data at kt ==!
CALL dta_tsd( kt, 'dmp', zts_dta ) ! read and interpolates T-S data at kt
!
SELECT CASE ( nn_zdmp ) !== type of damping ==!
!
CASE( 0 ) !* newtonian damping throughout the water column *!
DO jn = 1, jpts
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs) &
& + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) )
END_3D
END DO
!
CASE ( 1 ) !* no damping in the turbocline (avt > 5 cm2/s) *!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
IF( avt(ji,jj,jk) <= avt_c ) THEN
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) &
& + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) &
& + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
ENDIF
END_3D
!
CASE ( 2 ) !* no damping in the mixed layer *!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) &
& + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) &
& + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
ENDIF
END_3D
!
END SELECT
!
! outputs (clem trunk)
IF( iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN
ALLOCATE( zwrk(A2D(nn_hls),jpk) ) ! Needed to handle expressions containing e3t when using key_qco or key_linssh
zwrk(:,:,:) = 0._wp
IF( iom_use('hflx_dmp_cea') ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
zwrk(ji,jj,jk) = ( pts(ji,jj,jk,jp_tem,Krhs) - ztrdts(ji,jj,jk,jp_tem) ) * e3t(ji,jj,jk,Kmm)
END_3D
CALL iom_put('hflx_dmp_cea', SUM( zwrk(:,:,:), dim=3 ) * rcp * rho0 ) ! W/m2
ENDIF
IF( iom_use('sflx_dmp_cea') ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
zwrk(ji,jj,jk) = ( pts(ji,jj,jk,jp_sal,Krhs) - ztrdts(ji,jj,jk,jp_sal) ) * e3t(ji,jj,jk,Kmm)
END_3D
CALL iom_put('sflx_dmp_cea', SUM( zwrk(:,:,:), dim=3 ) * rho0 ) ! g/m2/s
ENDIF
DEALLOCATE( zwrk )
ENDIF
!
IF( l_trdtra ) THEN ! trend diagnostic
ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) - ztrdts(:,:,:,:)
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
DEALLOCATE( ztrdts )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp - Ta: ', mask1=tmask, &
& tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
!
IF( ln_timing ) CALL timing_stop('tra_dmp')
!
END SUBROUTINE tra_dmp
SUBROUTINE tra_dmp_init
!!----------------------------------------------------------------------
!! *** ROUTINE tra_dmp_init ***
!!
!! ** Purpose : Initialization for the newtonian damping
!!
!! ** Method : read the namtra_dmp namelist and check the parameters
!!----------------------------------------------------------------------
INTEGER :: ios, imask ! local integers
!
NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto
!!----------------------------------------------------------------------
!
READ ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in reference namelist' )
!
READ ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist' )
IF(lwm) WRITE ( numond, namtra_dmp )
!
IF(lwp) THEN ! Namelist print
WRITE(numout,*)
WRITE(numout,*) 'tra_dmp_init : T and S newtonian relaxation'
WRITE(numout,*) '~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namtra_dmp : set relaxation parameters'
WRITE(numout,*) ' Apply relaxation or not ln_tradmp = ', ln_tradmp
WRITE(numout,*) ' mixed layer damping option nn_zdmp = ', nn_zdmp
WRITE(numout,*) ' Damping file name cn_resto = ', cn_resto
WRITE(numout,*)
ENDIF
!
IF( ln_tradmp ) THEN
! ! Allocate arrays
IF( tra_dmp_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' )
!
SELECT CASE (nn_zdmp) ! Check values of nn_zdmp
CASE ( 0 ) ; IF(lwp) WRITE(numout,*) ' tracer damping as specified by mask'
CASE ( 1 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the mixing layer (kz > 5 cm2/s)'
CASE ( 2 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the mixed layer'
CASE DEFAULT
CALL ctl_stop('tra_dmp_init : wrong value of nn_zdmp')
END SELECT
!
!!TG: Initialisation of dtatsd - Would it be better to have dmpdta routine
! so can damp to something other than intitial conditions files?
!!gm: In principle yes. Nevertheless, we can't anticipate demands that have never been formulated.
IF( .NOT.ln_tsd_dmp ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout, *) ' read T-S data not initialized, we force ln_tsd_dmp=T'
CALL dta_tsd_init( ld_tradmp=ln_tradmp ) ! forces the initialisation of T-S data
ENDIF
! ! Read in mask from file
CALL iom_open ( cn_resto, imask)
CALL iom_get ( imask, jpdom_auto, 'resto', resto )
CALL iom_close( imask )
ENDIF
!
END SUBROUTINE tra_dmp_init
!!======================================================================
END MODULE tradmp
tests/ISOMIP/EXPREF/ISOMIP_mlt.png deleted 100644 → 0
View file @ 0b9e2937
tests/ISOMIP/EXPREF/ISOMIP_mlt.png

217 KiB

tests/ISOMIP/EXPREF/ISOMIP_moc.png deleted 100644 → 0
View file @ 0b9e2937
tests/ISOMIP/EXPREF/ISOMIP_moc.png

496 KiB