Changeset 15128

Timestamp:

2021-07-18T16:30:44+02:00 (3 years ago)

Author:

dbruciaferri

Message:

add 3D computation of model levels slope parameter

Location:

NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src

Files:

• domain.f90 (modified) (2 diffs)
• domwri.f90 (modified) (2 diffs)

NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domain.f90

@@ -385 +385 @@
       !!                              and depths (ln_e3_dep=F) 
       !!----------------------------------------------------------------------
-      INTEGER           ::   ji, jj, jk   ! dummy loop indices
-      INTEGER           ::   izco, izps, isco, icav
-      INTEGER           ::   inum     ! temprary units for 'domain_cfg.nc' file
-      CHARACTER(len=21) ::   clnam    ! filename (mesh and mask informations)
-      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! workspace
+      INTEGER                          ::   ji, jj, jk             ! dummy loop indices
+      INTEGER                          ::   izco, izps, isco, icav
+      INTEGER                          ::   inum                   ! temporary units for 
+                                                                   ! 'domain_cfg.nc' file
+      CHARACTER(len=21)                ::   clnam                  ! filename (mesh and 
+                                                                   ! mask informations)
+      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                    ! workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d                    ! workspace
       !!----------------------------------------------------------------------
       !
@@ -493 +496 @@
          CALL dom_stiff( z2d )
          CALL iom_rstput( 0, 0, inum, 'stiffness', z2d )        !    ! Max. grid stiffness ratio
+         CALL dom_stiff_3D( z3d )
+         CALL iom_rstput( 0, 0, inum, 'stiff3D', z3d )  ! 3D stiffness ratio
       ENDIF
       !

NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domwri.f90

@@ -30 +30 @@
    PUBLIC   dom_wri_coordinate   ! routine called by domhgr.F90
    PUBLIC   dom_stiff            ! routine called by inidom.F90
+   PUBLIC   dom_stiff_3d         ! routine called by domain.F90
 
    !!----------------------------------------------------------------------
@@ -480 +481 @@
    END SUBROUTINE dom_stiff
 
+
+   SUBROUTINE dom_stiff_3D( px1 )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_stiff  ***
+      !!                     
+      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
+      !!
+      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
+      !!                Save the maximum in the vertical direction
+      !!                (this number is only relevant in s-coordinates)
+      !!
+      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
+      !
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zrxmax
+      REAL(wp), DIMENSION(4) ::   zr1
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zx1
+      !!----------------------------------------------------------------------
+      zx1(:,:,:) = 0._wp
+      zrxmax     = 0._wp
+      zr1(:)     = 0._wp
+      !
+      DO ji = 2, jpim1
+         DO jj = 2, jpjm1
+            DO jk = 1, jpkm1
+!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation ....
+!!             especially since it is gde3w which is used to compute the pressure gradient
+!!             furthermore, I think gdept_0 should be used below instead of w-point in the numerator
+!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
+               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               &
+                    &           +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )            &
+                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
+                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
+               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
+                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
+                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
+                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
+               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
+                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
+                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
+                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
+               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
+                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
+                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
+                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
+               zrxmax = MAXVAL( zr1(1:4) )
+               zx1(ji,jj,jk) = zrxmax
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( zx1, 'T', 1. )
+      !
+      IF( PRESENT( px1 ) )    px1 = zx1
+      !
+   END SUBROUTINE dom_stiff_3D
+
    !!======================================================================
 END MODULE domwri