New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 3618 for branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 – NEMO

Ignore:
Timestamp:
2012-11-20T19:20:57+01:00 (10 years ago)
Author:
johnsiddorn
Message:

updates to reflect reviewers comments on coding standards and documentation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r3600 r3618  
    1515   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option 
    1616   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level 
    17    !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn adn Furner stretching functio 
     17   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function 
    1818   !!---------------------------------------------------------------------- 
    1919 
     
    2828   !!       zgr_zps      : z-coordinate with partial steps 
    2929   !!       zgr_sco      : s-coordinate 
    30    !!       fssig        : sigma coordinate non-dimensional function 
    31    !!       fgamma       : Siddorn and Furner stretching function 
    32    !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!) 
     30   !!       fssig        : tanh stretch function 
     31   !!       fssig1       : Song and Haidvogel 1994 stretch function 
     32   !!       fgamma       : Siddorn and Furner 2012 stretching function 
    3333   !!--------------------------------------------------------------------- 
    3434   USE oce               ! ocean variables 
     
    5151   LOGICAL  ::   ln_s_sh94   = .false.      ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) 
    5252   LOGICAL  ::   ln_s_sf12   = .true.       ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T) 
    53    LOGICAL  ::   ln_sigcrit  = .false.      ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch  
    5453   ! 
    5554   REAL(wp) ::   rn_sbot_min =  300._wp     ! minimum depth of s-bottom surface (>0) (m) 
     
    6362   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 
    6463   ! Siddorn and Furner stretching parameters 
     64   LOGICAL  ::   ln_sigcrit  = .false.      ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch  
    6565   REAL(wp) ::   rn_alpha    =    4.4_wp    ! control parameter ( > 1 stretch towards surface, < 1 towards seabed) 
    6666   REAL(wp) ::   rn_efold    =    0.0_wp    !  efold length scale for transition to stretched coord 
     
    10651065      !!            hbatv = mj( hbatt ) 
    10661066      !!            hbatf = mi( mj( hbatt ) ) 
    1067       !!          - Compute gsigt, gsigw, esigt, esigw from an analytical 
     1067      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical 
    10681068      !!         function and its derivative given as function. 
    1069       !!            gsigt(k) = fssig (k    ) 
    1070       !!            gsigw(k) = fssig (k-0.5) 
    1071       !!            esigt(k) = fsdsig(k    ) 
    1072       !!            esigw(k) = fsdsig(k-0.5) 
     1069      !!            z_gsigt(k) = fssig (k    ) 
     1070      !!            z_gsigw(k) = fssig (k-0.5) 
     1071      !!            z_esigt(k) = fsdsig(k    ) 
     1072      !!            z_esigw(k) = fsdsig(k-0.5) 
    10731073      !!      Three options for stretching are give, and they can be modified 
    10741074      !!      following the users requirements. Nevertheless, the output as 
     
    11251125         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb 
    11261126         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12 
     1127         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit 
    11271128         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients' 
    11281129         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha 
     
    13831384         &                                                       ' MAX ', MAXVAL( mbathy(:,:) ) 
    13841385 
    1385       !                                               ! ============= 
    1386       IF(lwp) THEN                                    ! Control print 
    1387          !                                            ! ============= 
    1388          WRITE(numout,*)  
    1389          WRITE(numout,*) ' domzgr: vertical coefficients for model level' 
    1390          WRITE(numout, "(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')" ) 
    1391          WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk ) 
    1392       ENDIF 
    13931386      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain 
    13941387         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)   ), ' MAX ', MAXVAL( mbathy(:,:) ) 
     
    15081501      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
    15091502      ! 
    1510       REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 
    1511       REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3            
    1512  
    1513       CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
    1514       CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
    1515  
    1516       gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp 
    1517       esigt3  = 0._wp   ;   esigw3  = 0._wp  
    1518       esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp 
    1519       esigwu3 = 0._wp   ;   esigwv3 = 0._wp 
     1503      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     1504      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3            
     1505 
     1506      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      ) 
     1507      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
     1508 
     1509      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp 
     1510      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp  
     1511      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp 
     1512      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp 
    15201513 
    15211514      DO ji = 1, jpi 
     
    15241517            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
    15251518               DO jk = 1, jpk 
    1526                   gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 
    1527                   gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb ) 
     1519                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 
     1520                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb ) 
    15281521               END DO 
    15291522            ELSE ! shallow water, uniform sigma 
    15301523               DO jk = 1, jpk 
    1531                   gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp) 
    1532                   gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 
     1524                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp) 
     1525                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 
    15331526                  END DO 
    15341527            ENDIF 
    15351528            ! 
    15361529            DO jk = 1, jpkm1 
    1537                esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 
    1538                esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 
    1539             END DO 
    1540             esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) ) 
    1541             esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 
     1530               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 
     1531               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 
     1532            END DO 
     1533            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) ) 
     1534            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) 
    15421535            ! 
    15431536            ! Coefficients for vertical depth as the sum of e3w scale factors 
    1544             gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 
     1537            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1) 
    15451538            DO jk = 2, jpk 
    1546                gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
     1539               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 
    15471540            END DO 
    15481541            ! 
     
    15501543               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
    15511544               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
    1552                gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
    1553                gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
    1554                gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 
     1545               gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
     1546               gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
     1547               gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 
    15551548            END DO 
    15561549           ! 
     
    15611554         DO jj = 1, jpjm1 
    15621555            DO jk = 1, jpk 
    1563                esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   & 
     1556               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   & 
    15641557                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1565                esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   & 
     1558               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   & 
    15661559                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1567                esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     & 
    1568                   &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   & 
     1560               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     & 
     1561                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   & 
    15691562                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    1570                esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   & 
     1563               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   & 
    15711564                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1572                esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   & 
     1565               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   & 
    15731566                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    15741567               ! 
    1575                e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
    1576                e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
    1577                e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
    1578                e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1568               e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1569               e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1570               e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1571               e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
    15791572               ! 
    1580                e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
    1581                e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
    1582                e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1573               e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1574               e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1575               e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
    15831576            END DO 
    15841577        END DO 
    15851578      END DO 
    15861579 
    1587       CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
    1588       CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
     1580      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      ) 
     1581      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
    15891582 
    15901583   END SUBROUTINE s_sh94 
     
    16091602      ! 
    16101603      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
    1611       REAL(wp) ::   fsmth          ! smoothing around critical depth 
    1612       REAL(wp) ::   zss, zbb       ! Surface and bottom cell thickness in sigma space 
    1613       ! 
    1614       REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 
    1615       REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3            
    1616  
    1617       ! 
    1618       CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
    1619       CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
    1620  
    1621       gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp 
    1622       esigt3  = 0._wp   ;   esigw3  = 0._wp  
    1623       esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp 
    1624       esigwu3 = 0._wp   ;   esigwv3 = 0._wp 
     1604      REAL(wp) ::   zsmth               ! smoothing around critical depth 
     1605      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space 
     1606      ! 
     1607      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     1608      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3            
     1609 
     1610      ! 
     1611      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      ) 
     1612      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
     1613 
     1614      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp 
     1615      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp  
     1616      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp 
     1617      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp 
    16251618 
    16261619      DO ji = 1, jpi 
     
    16291622          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 
    16301623               
    1631               zbb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,. 
     1624              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,. 
    16321625                                                     ! could be changed by users but care must be taken to do so carefully 
    1633               zbb = 1.0_wp-(zbb/hbatt(ji,jj)) 
     1626              zzb = 1.0_wp-(zzb/hbatt(ji,jj)) 
    16341627             
    1635               zss = rn_zs / hbatt(ji,jj)  
     1628              zzs = rn_zs / hbatt(ji,jj)  
    16361629               
    16371630              IF (rn_efold /= 0.0_wp) THEN 
    1638                 fsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 
     1631                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 
    16391632              ELSE 
    1640                 fsmth = 1.0_wp  
     1633                zsmth = 1.0_wp  
    16411634              ENDIF 
    16421635                
    16431636              DO jk = 1, jpk 
    1644                 gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp) 
    1645                 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp) 
     1637                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp) 
     1638                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp) 
    16461639              ENDDO 
    1647               gsigw3(ji,jj,:) = fgamma( gsigw3(ji,jj,:), zbb, zss, fsmth  ) 
    1648               gsigt3(ji,jj,:) = fgamma( gsigt3(ji,jj,:), zbb, zss, fsmth  ) 
     1640              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  ) 
     1641              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  ) 
    16491642  
    16501643          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma 
    16511644 
    16521645            DO jk = 1, jpk 
    1653               gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp) 
    1654               gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 
     1646              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp) 
     1647              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 
    16551648            END DO 
    16561649 
     
    16581651 
    16591652            DO jk = 1, jpk 
    1660               gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
    1661               gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
     1653              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
     1654              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
    16621655            END DO 
    16631656 
     
    16651658 
    16661659          DO jk = 1, jpkm1 
    1667              esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 
    1668              esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 
     1660             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 
     1661             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 
    16691662          END DO 
    1670           esigw3(ji,jj,1  ) = 2.0_wp * (gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  )) 
    1671           esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk)) 
     1663          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  )) 
     1664          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) 
    16721665 
    16731666          ! Coefficients for vertical depth as the sum of e3w scale factors 
    1674           gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1) 
     1667          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 
    16751668          DO jk = 2, jpk 
    1676              gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
     1669             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 
    16771670          END DO 
    16781671 
    16791672          DO jk = 1, jpk 
    1680              gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigt3(ji,jj,jk) 
    1681              gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigw3(ji,jj,jk) 
    1682              gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsi3w3(ji,jj,jk) 
     1673             gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 
     1674             gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 
     1675             gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 
    16831676          END DO 
    16841677 
     
    16901683 
    16911684          DO jk = 1, jpk 
    1692                 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / & 
     1685                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 
    16931686                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1694                 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / & 
     1687                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 
    16951688                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1696                 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) +  & 
    1697                                       hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / & 
     1689                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  & 
     1690                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 
    16981691                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    1699                 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / & 
     1692                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 
    17001693                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1701                 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / & 
     1694                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 
    17021695                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    17031696 
    1704              e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*esigt3(ji,jj,jk) 
    1705              e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*esigtu3(ji,jj,jk) 
    1706              e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*esigtv3(ji,jj,jk) 
    1707              e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*esigtf3(ji,jj,jk) 
     1697             e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 
     1698             e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk) 
     1699             e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk) 
     1700             e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 
    17081701             ! 
    1709              e3w(ji,jj,jk)=hbatt(ji,jj)*esigw3(ji,jj,jk) 
    1710              e3uw(ji,jj,jk)=hbatu(ji,jj)*esigwu3(ji,jj,jk) 
    1711              e3vw(ji,jj,jk)=hbatv(ji,jj)*esigwv3(ji,jj,jk) 
     1702             e3w(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 
     1703             e3uw(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 
     1704             e3vw(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) 
    17121705          END DO 
    17131706 
    17141707        ENDDO 
    17151708      ENDDO 
    1716  
    1717       CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
    1718       CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
     1709      !                                               ! ============= 
     1710 
     1711      CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      ) 
     1712      CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
    17191713 
    17201714   END SUBROUTINE s_sf12 
     
    17351729      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
    17361730 
     1731      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 
     1732      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 
     1733 
     1734      CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      ) 
     1735      CALL wrk_alloc( jpk, z_esigt, z_esigw                                               ) 
     1736 
     1737      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp 
     1738      z_esigt  = 0._wp   ;   z_esigw  = 0._wp  
     1739 
    17371740      DO jk = 1, jpk 
    1738         gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 
    1739         gsigt(jk) = -fssig( REAL(jk,wp)        ) 
    1740       END DO 
    1741       IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk) 
     1741        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 
     1742        z_gsigt(jk) = -fssig( REAL(jk,wp)        ) 
     1743      END DO 
     1744      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk) 
    17421745      ! 
    17431746      ! Coefficients for vertical scale factors at w-, t- levels 
     
    17451748!!gm        or betteroffer the 2 possibilities.... 
    17461749      DO jk = 1, jpkm1 
    1747          esigt(jk  ) = gsigw(jk+1) - gsigw(jk) 
    1748          esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 
    1749       END DO 
    1750       esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) )  
    1751       esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 
     1750         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk) 
     1751         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk) 
     1752      END DO 
     1753      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) )  
     1754      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) ) 
    17521755      ! 
    17531756      ! Coefficients for vertical depth as the sum of e3w scale factors 
    1754       gsi3w(1) = 0.5_wp * esigw(1) 
     1757      z_gsi3w(1) = 0.5_wp * z_esigw(1) 
    17551758      DO jk = 2, jpk 
    1756          gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 
     1759         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk) 
    17571760      END DO 
    17581761!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
     
    17601763         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
    17611764         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
    1762          gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 
    1763          gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 
    1764          gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 
     1765         gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 
     1766         gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 
     1767         gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 
    17651768      END DO 
    17661769!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
     
    17681771         DO ji = 1, jpi 
    17691772            DO jk = 1, jpk 
    1770               e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
    1771               e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
    1772               e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
    1773               e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 
     1773              e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
     1774              e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
     1775              e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
     1776              e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 
    17741777              ! 
    1775               e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
    1776               e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
    1777               e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
    1778             END DO 
    1779          END DO 
    1780       END DO 
     1778              e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
     1779              e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
     1780              e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
     1781            END DO 
     1782         END DO 
     1783      END DO 
     1784 
     1785      CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w                                      ) 
     1786      CALL wrk_dealloc( jpk, z_esigt, z_esigw                                               ) 
     1787 
    17811788   END SUBROUTINE s_tanh 
    17821789 
     
    18321839 
    18331840 
    1834    FUNCTION fgamma( pk1, Zbb, Zss, F ) RESULT( gam ) 
     1841   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma ) 
    18351842      !!---------------------------------------------------------------------- 
    18361843      !!                 ***  ROUTINE fgamma  *** 
     
    18491856      !! Reference  :   Siddorn and Furner, in prep 
    18501857      !!---------------------------------------------------------------------- 
    1851       REAL(wp), INTENT(in   ) ::   pk1(jpk)   ! continuous "k" coordinate 
    1852       REAL(wp)                ::   gam(jpk)   ! stretched coordinate 
    1853       REAL(wp), INTENT(in   ) ::   Zbb      ! Bottom box depth 
    1854       REAL(wp), INTENT(in   ) ::   Zss      ! surface box depth 
    1855       REAL(wp), INTENT(in   ) ::   F        ! Smoothing parameter 
    1856       REAL(wp)                ::   a1,a2,a3 ! local variables 
    1857       REAL(wp)                ::   n1,n2    ! local variables 
    1858       REAL(wp)                ::   A,B,X    ! local variables 
     1858      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate 
     1859      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate 
     1860      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth 
     1861      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth 
     1862      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter 
     1863      REAL(wp)                ::   za1,za2,za3    ! local variables 
     1864      REAL(wp)                ::   zn1,zn2        ! local variables 
     1865      REAL(wp)                ::   za,zb,zx       ! local variables 
    18591866      integer                 ::   jk 
    18601867      !!---------------------------------------------------------------------- 
    18611868      ! 
    18621869 
    1863       n1  =  1./(jpk-1.) 
    1864       n2  =  1. -  n1 
    1865  
    1866       a1 = (rn_alpha+2.0_wp)*n1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n1**(rn_alpha+2.0_wp)  
    1867       a2 = (rn_alpha+2.0_wp)*n2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n2**(rn_alpha+2.0_wp) 
    1868       a3 = ( n2**3.0_wp - a2)/( n1**3.0_wp - a1) 
     1870      zn1  =  1./(jpk-1.) 
     1871      zn2  =  1. -  zn1 
     1872 
     1873      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp)  
     1874      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 
     1875      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 
    18691876      
    1870       A = Zbb - a3*(Zss-a1)-a2 
    1871       A = A/( n2-0.5_wp*(a2+n2**2.0_wp) - a3*(n1-0.5_wp*(a1+n1**2.0_wp) ) ) 
    1872       B = (Zss - a1 - A*( n1-0.5_wp*(a1+n1**2.0_wp ) ) ) / (n1**3.0_wp - a1) 
    1873       X = 1.0_wp-A/2.0_wp-B 
     1877      za = pzb - za3*(pzs-za1)-za2 
     1878      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 
     1879      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 
     1880      zx = 1.0_wp-za/2.0_wp-zb 
    18741881  
    18751882      DO jk = 1, jpk 
    1876         gam(jk) = A*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+B*pk1(jk)**3.0_wp + X*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 
    1877         gam(jk) = gam(jk)*F+pk1(jk)*(1.0_wp-F) 
     1883        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 
     1884        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 
    18781885      ENDDO  
    18791886 
Note: See TracChangeset for help on using the changeset viewer.