Ignore:
Timestamp:
2012-11-27T15:42:24+01:00 (9 years ago)
Author:
rblod
Message:

First commit of the final branch for 2012 (future nemo_3_5), see ticket #1028

File:
1 edited

Legend:

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

    r3632 r3680  
    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 and Furner stretching function 
    1718   !!---------------------------------------------------------------------- 
    1819 
     
    2728   !!       zgr_zps      : z-coordinate with partial steps 
    2829   !!       zgr_sco      : s-coordinate 
    29    !!       fssig        : sigma coordinate non-dimensional function 
    30    !!       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 
    3133   !!--------------------------------------------------------------------- 
    3234   USE oce               ! ocean variables 
     
    4749 
    4850   !                                       !!* Namelist namzgr_sco * 
     51   LOGICAL  ::   ln_s_sh94   = .false.      ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) 
     52   LOGICAL  ::   ln_s_sf12   = .true.       ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T) 
     53   ! 
    4954   REAL(wp) ::   rn_sbot_min =  300._wp     ! minimum depth of s-bottom surface (>0) (m) 
    5055   REAL(wp) ::   rn_sbot_max = 5250._wp     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     56   REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1) 
     57   REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for transition from sigma to stretched coordinates 
     58   ! Song and Haidvogel 1994 stretching parameters 
    5159   REAL(wp) ::   rn_theta    =    6.00_wp   ! surface control parameter (0<=rn_theta<=20) 
    5260   REAL(wp) ::   rn_thetb    =    0.75_wp   ! bottom control parameter  (0<=rn_thetb<= 1) 
    53    REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1) 
    54    LOGICAL  ::   ln_s_sigma  = .false.      ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T) 
    55    REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter for song and haidvogel stretching 
     61   REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter  
    5662   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 
    57    REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for s-sigma coordinates 
     63   ! 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  
     65   REAL(wp) ::   rn_alpha    =    4.4_wp    ! control parameter ( > 1 stretch towards surface, < 1 towards seabed) 
     66   REAL(wp) ::   rn_efold    =    0.0_wp    !  efold length scale for transition to stretched coord 
     67   REAL(wp) ::   rn_zs       =    1.0_wp    !  depth of surface grid box 
     68                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 
     69   REAL(wp) ::   rn_zb_a     =    0.024_wp  !  bathymetry scaling factor for calculating Zb 
     70   REAL(wp) ::   rn_zb_b     =   -0.2_wp    !  offset for calculating Zb 
    5871 
    5972  !! * Substitutions 
     
    10341047   END SUBROUTINE zgr_zps 
    10351048 
    1036  
    1037    FUNCTION fssig( pk ) RESULT( pf ) 
    1038       !!---------------------------------------------------------------------- 
    1039       !!                 ***  ROUTINE eos_init  *** 
    1040       !!        
    1041       !! ** Purpose :   provide the analytical function in s-coordinate 
    1042       !!           
    1043       !! ** Method  :   the function provide the non-dimensional position of 
    1044       !!                T and W (i.e. between 0 and 1) 
    1045       !!                T-points at integer values (between 1 and jpk) 
    1046       !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
    1047       !!---------------------------------------------------------------------- 
    1048       REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate 
    1049       REAL(wp)             ::   pf   ! sigma value 
    1050       !!---------------------------------------------------------------------- 
    1051       ! 
    1052       pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   & 
    1053          &     - TANH( rn_thetb * rn_theta                                )  )   & 
    1054          & * (   COSH( rn_theta                           )                      & 
    1055          &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              & 
    1056          & / ( 2._wp * SINH( rn_theta ) ) 
    1057       ! 
    1058    END FUNCTION fssig 
    1059  
    1060  
    1061    FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 
    1062       !!---------------------------------------------------------------------- 
    1063       !!                 ***  ROUTINE eos_init  *** 
    1064       !! 
    1065       !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate 
    1066       !! 
    1067       !! ** Method  :   the function provides the non-dimensional position of 
    1068       !!                T and W (i.e. between 0 and 1) 
    1069       !!                T-points at integer values (between 1 and jpk) 
    1070       !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
    1071       !!---------------------------------------------------------------------- 
    1072       REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate 
    1073       REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient 
    1074       REAL(wp)             ::   pf1   ! sigma value 
    1075       !!---------------------------------------------------------------------- 
    1076       ! 
    1077       IF ( rn_theta == 0 ) then      ! uniform sigma 
    1078          pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 
    1079       ELSE                        ! stretched sigma 
    1080          pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              & 
    1081             &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  & 
    1082             &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  ) 
    1083       ENDIF 
    1084       ! 
    1085    END FUNCTION fssig1 
    1086  
    1087  
    10881049   SUBROUTINE zgr_sco 
    10891050      !!---------------------------------------------------------------------- 
     
    11041065      !!            hbatv = mj( hbatt ) 
    11051066      !!            hbatf = mi( mj( hbatt ) ) 
    1106       !!          - Compute gsigt, gsigw, esigt, esigw from an analytical 
     1067      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical 
    11071068      !!         function and its derivative given as function. 
    1108       !!            gsigt(k) = fssig (k    ) 
    1109       !!            gsigw(k) = fssig (k-0.5) 
    1110       !!            esigt(k) = fsdsig(k    ) 
    1111       !!            esigw(k) = fsdsig(k-0.5) 
    1112       !!      This routine is given as an example, it must be modified 
    1113       !!      following the user s desiderata. nevertheless, the output as 
     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) 
     1073      !!      Three options for stretching are give, and they can be modified 
     1074      !!      following the users requirements. Nevertheless, the output as 
    11141075      !!      well as the way to compute the model levels and scale factors 
    1115       !!      must be respected in order to insure second order a!!uracy 
     1076      !!      must be respected in order to insure second order accuracy 
    11161077      !!      schemes. 
    11171078      !! 
    1118       !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 
     1079      !!      The three methods for stretching available are: 
     1080      !!  
     1081      !!           s_sh94 (Song and Haidvogel 1994) 
     1082      !!                a sinh/tanh function that allows sigma and stretched sigma 
     1083      !! 
     1084      !!           s_sf12 (Siddorn and Furner 2012?) 
     1085      !!                allows the maintenance of fixed surface and or 
     1086      !!                bottom cell resolutions (cf. geopotential coordinates)  
     1087      !!                within an analytically derived stretched S-coordinate framework. 
     1088      !!  
     1089      !!          s_tanh  (Madec et al 1996) 
     1090      !!                a cosh/tanh function that gives stretched coordinates         
     1091      !! 
    11191092      !!---------------------------------------------------------------------- 
    11201093      ! 
    11211094      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument 
    11221095      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers 
    1123       REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper   ! temporary scalars 
     1096      REAL(wp) ::   zrmax, ztaper   ! temporary scalars 
    11241097      ! 
    11251098      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 
    1126       REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 
    1127       REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3            
    1128  
    1129       NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc 
    1130       !!---------------------------------------------------------------------- 
     1099 
     1100      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 
     1101                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 
     1102     !!---------------------------------------------------------------------- 
    11311103      ! 
    11321104      IF( nn_timing == 1 )  CALL timing_start('zgr_sco') 
    11331105      ! 
    11341106      CALL wrk_alloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           ) 
    1135       CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
    1136       CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
    11371107      ! 
    11381108      REWIND( numnam )                       ! Read Namelist namzgr_sco : sigma-stretching parameters 
     
    11441114         WRITE(numout,*) '~~~~~~~~~~~' 
    11451115         WRITE(numout,*) '   Namelist namzgr_sco' 
    1146          WRITE(numout,*) '      sigma-stretching coeffs ' 
    1147          WRITE(numout,*) '      maximum depth of s-bottom surface (>0)       rn_sbot_max   = ' ,rn_sbot_max 
    1148          WRITE(numout,*) '      minimum depth of s-bottom surface (>0)       rn_sbot_min   = ' ,rn_sbot_min 
    1149          WRITE(numout,*) '      surface control parameter (0<=rn_theta<=20)  rn_theta      = ', rn_theta 
    1150          WRITE(numout,*) '      bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ', rn_thetb 
    1151          WRITE(numout,*) '      maximum cut-off r-value allowed              rn_rmax       = ', rn_rmax 
    1152          WRITE(numout,*) '      Hybrid s-sigma-coordinate                    ln_s_sigma    = ', ln_s_sigma 
    1153          WRITE(numout,*) '      stretching parameter (song and haidvogel)    rn_bb         = ', rn_bb 
    1154          WRITE(numout,*) '      Critical depth                               rn_hc         = ', rn_hc 
    1155       ENDIF 
    1156  
    1157       gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp 
    1158       esigt3  = 0._wp   ;   esigw3  = 0._wp  
    1159       esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp 
    1160       esigwu3 = 0._wp   ;   esigwv3 = 0._wp 
     1116         WRITE(numout,*) '     stretching coeffs ' 
     1117         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max 
     1118         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min 
     1119         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc 
     1120         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax 
     1121         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94 
     1122         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients' 
     1123         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta 
     1124         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb 
     1125         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb 
     1126         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 
     1128         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients' 
     1129         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha 
     1130         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold 
     1131         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs 
     1132         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a 
     1133         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b 
     1134         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 
     1135      ENDIF 
    11611136 
    11621137      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate 
     
    13521327      ! non-dimensional "sigma" for model level depth at w- and t-levels 
    13531328 
    1354       IF( ln_s_sigma ) THEN        ! Song and Haidvogel style stretched sigma for depths 
    1355          !                         ! below rn_hc, with uniform sigma in shallower waters 
    1356          DO ji = 1, jpi 
    1357             DO jj = 1, jpj 
    1358  
    1359                IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
    1360                   DO jk = 1, jpk 
    1361                      gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 
    1362                      gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb ) 
    1363                   END DO 
    1364                ELSE ! shallow water, uniform sigma 
    1365                   DO jk = 1, jpk 
    1366                      gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp) 
    1367                      gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 
    1368                   END DO 
    1369                ENDIF 
    1370                IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw3 1 jpk    ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 
    1371                ! 
    1372                DO jk = 1, jpkm1 
    1373                   esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 
    1374                   esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 
    1375                END DO 
    1376                esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) ) 
    1377                esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 
    1378                ! 
    1379                ! Coefficients for vertical depth as the sum of e3w scale factors 
    1380                gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 
    1381                DO jk = 2, jpk 
    1382                   gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
    1383                END DO 
    1384                ! 
    1385                DO jk = 1, jpk 
    1386                   zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
    1387                   zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
    1388                   gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
    1389                   gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
    1390                   gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 
    1391                END DO 
    1392                ! 
    1393             END DO   ! for all jj's 
    1394          END DO    ! for all ji's 
    1395  
    1396          DO ji = 1, jpim1 
    1397             DO jj = 1, jpjm1 
    1398                DO jk = 1, jpk 
    1399                   esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   & 
    1400                      &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1401                   esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   & 
    1402                      &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1403                   esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     & 
    1404                      &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   & 
    1405                      &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    1406                   esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   & 
    1407                      &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1408                   esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   & 
    1409                      &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1410                   ! 
    1411                   e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1412                   e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1413                   e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1414                   e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1415                   ! 
    1416                   e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1417                   e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1418                   e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1419                END DO 
    1420             END DO 
    1421          END DO 
    1422  
    1423          CALL lbc_lnk( e3t , 'T', 1._wp ) 
    1424          CALL lbc_lnk( e3u , 'U', 1._wp ) 
    1425          CALL lbc_lnk( e3v , 'V', 1._wp ) 
    1426          CALL lbc_lnk( e3f , 'F', 1._wp ) 
    1427          CALL lbc_lnk( e3w , 'W', 1._wp ) 
    1428          CALL lbc_lnk( e3uw, 'U', 1._wp ) 
    1429          CALL lbc_lnk( e3vw, 'V', 1._wp ) 
    1430  
    1431          ! 
    1432       ELSE   ! not ln_s_sigma 
    1433          ! 
    1434          DO jk = 1, jpk 
    1435            gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 
    1436            gsigt(jk) = -fssig( REAL(jk,wp)        ) 
    1437          END DO 
    1438          IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk) 
    1439          ! 
    1440          ! Coefficients for vertical scale factors at w-, t- levels 
    1441 !!gm bug :  define it from analytical function, not like juste bellow.... 
    1442 !!gm        or betteroffer the 2 possibilities.... 
    1443          DO jk = 1, jpkm1 
    1444             esigt(jk  ) = gsigw(jk+1) - gsigw(jk) 
    1445             esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 
    1446          END DO 
    1447          esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) )  
    1448          esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 
    1449  
    1450 !!gm  original form 
    1451 !!org DO jk = 1, jpk 
    1452 !!org    esigt(jk)=fsdsig( FLOAT(jk)     ) 
    1453 !!org    esigw(jk)=fsdsig( FLOAT(jk)-0.5 ) 
    1454 !!org END DO 
    1455 !!gm 
    1456          ! 
    1457          ! Coefficients for vertical depth as the sum of e3w scale factors 
    1458          gsi3w(1) = 0.5_wp * esigw(1) 
    1459          DO jk = 2, jpk 
    1460             gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 
    1461          END DO 
    1462 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
    1463          DO jk = 1, jpk 
    1464             zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
    1465             zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
    1466             gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 
    1467             gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 
    1468             gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 
    1469          END DO 
    1470 !!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
    1471          DO jj = 1, jpj 
    1472             DO ji = 1, jpi 
    1473                DO jk = 1, jpk 
    1474                  e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
    1475                  e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
    1476                  e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
    1477                  e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 
    1478                  ! 
    1479                  e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
    1480                  e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
    1481                  e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
    1482                END DO 
    1483             END DO 
    1484          END DO 
    1485          ! 
    1486       ENDIF ! ln_s_sigma 
    1487  
    1488  
     1329 
     1330!======================================================================== 
     1331! Song and Haidvogel  1994 (ln_s_sh94=T) 
     1332! Siddorn and Furner 2012 (ln_sf12=T) 
     1333! or  tanh function       (both false)                     
     1334!======================================================================== 
     1335      IF      ( ln_s_sh94 ) THEN  
     1336                           CALL s_sh94() 
     1337      ELSE IF ( ln_s_sf12 ) THEN 
     1338                           CALL s_sf12() 
     1339      ELSE                  
     1340                           CALL s_tanh() 
     1341      ENDIF  
     1342 
     1343      CALL lbc_lnk( e3t , 'T', 1._wp ) 
     1344      CALL lbc_lnk( e3u , 'U', 1._wp ) 
     1345      CALL lbc_lnk( e3v , 'V', 1._wp ) 
     1346      CALL lbc_lnk( e3f , 'F', 1._wp ) 
     1347      CALL lbc_lnk( e3w , 'W', 1._wp ) 
     1348      CALL lbc_lnk( e3uw, 'U', 1._wp ) 
     1349      CALL lbc_lnk( e3vw, 'V', 1._wp ) 
     1350 
     1351      fsdepw(:,:,:) = gdepw (:,:,:) 
     1352      fsde3w(:,:,:) = gdep3w(:,:,:) 
    14891353      ! 
    14901354      where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1.0 
     
    15201384         &                                                       ' MAX ', MAXVAL( mbathy(:,:) ) 
    15211385 
    1522       !                                               ! ============= 
    1523       IF(lwp) THEN                                    ! Control print 
    1524          !                                            ! ============= 
    1525          WRITE(numout,*)  
    1526          WRITE(numout,*) ' domzgr: vertical coefficients for model level' 
    1527          WRITE(numout, "(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')" ) 
    1528          WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk ) 
    1529       ENDIF 
    15301386      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain 
    15311387         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)   ), ' MAX ', MAXVAL( mbathy(:,:) ) 
     
    15441400            &                          ' w ', MAXVAL( fse3w (:,:,:) ) 
    15451401      ENDIF 
    1546       ! 
     1402      !  END DO 
    15471403      IF(lwp) THEN                                  ! selected vertical profiles 
    15481404         WRITE(numout,*) 
     
    15741430      ENDIF 
    15751431 
    1576 !!gm bug?  no more necessary?  if ! defined key_helsinki 
     1432!================================================================================ 
     1433! check the coordinate makes sense 
     1434!================================================================================ 
     1435      DO ji = 1, jpi 
     1436         DO jj = 1, jpj 
     1437 
     1438            IF( hbatt(ji,jj) > 0._wp) THEN 
     1439               DO jk = 1, mbathy(ji,jj) 
     1440                 ! check coordinate is monotonically increasing 
     1441                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
     1442                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1443                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1444                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 
     1445                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 
     1446                    CALL ctl_stop( ctmp1 ) 
     1447                 ENDIF 
     1448                 ! and check it has never gone negative 
     1449                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 
     1450                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
     1451                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1452                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
     1453                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
     1454                    CALL ctl_stop( ctmp1 ) 
     1455                 ENDIF 
     1456                 ! and check it never exceeds the total depth 
     1457                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1458                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1459                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1460                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
     1461                    CALL ctl_stop( ctmp1 ) 
     1462                 ENDIF 
     1463               END DO 
     1464 
     1465               DO jk = 1, mbathy(ji,jj)-1 
     1466                 ! and check it never exceeds the total depth 
     1467                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1468                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1469                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1470                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
     1471                    CALL ctl_stop( ctmp1 ) 
     1472                 ENDIF 
     1473               END DO 
     1474 
     1475            ENDIF 
     1476 
     1477         END DO 
     1478      END DO 
     1479      ! 
     1480      CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           ) 
     1481      ! 
     1482      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco') 
     1483      ! 
     1484   END SUBROUTINE zgr_sco 
     1485 
     1486!!====================================================================== 
     1487   SUBROUTINE s_sh94() 
     1488 
     1489      !!---------------------------------------------------------------------- 
     1490      !!                  ***  ROUTINE s_sh94  *** 
     1491      !!                      
     1492      !! ** Purpose :   stretch the s-coordinate system 
     1493      !! 
     1494      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994 
     1495      !!                mixed S/sigma coordinate 
     1496      !! 
     1497      !! Reference : Song and Haidvogel 1994.  
     1498      !!---------------------------------------------------------------------- 
     1499      ! 
     1500      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     1501      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     1502      ! 
     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 
     1513 
     1514      DO ji = 1, jpi 
     1515         DO jj = 1, jpj 
     1516 
     1517            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
     1518               DO jk = 1, jpk 
     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 ) 
     1521               END DO 
     1522            ELSE ! shallow water, uniform sigma 
     1523               DO jk = 1, jpk 
     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) 
     1526                  END DO 
     1527            ENDIF 
     1528            ! 
     1529            DO jk = 1, jpkm1 
     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) ) 
     1535            ! 
     1536            ! Coefficients for vertical depth as the sum of e3w scale factors 
     1537            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1) 
     1538            DO jk = 2, jpk 
     1539               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 
     1540            END DO 
     1541            ! 
     1542            DO jk = 1, jpk 
     1543               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
     1544               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
     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 ) 
     1548            END DO 
     1549           ! 
     1550         END DO   ! for all jj's 
     1551      END DO    ! for all ji's 
     1552 
     1553      DO ji = 1, jpim1 
     1554         DO jj = 1, jpjm1 
     1555            DO jk = 1, jpk 
     1556               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   & 
     1557                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1558               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   & 
     1559                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     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) )   & 
     1562                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     1563               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   & 
     1564                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1565               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   & 
     1566                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1567               ! 
     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) ) 
     1572               ! 
     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) ) 
     1576            END DO 
     1577        END DO 
     1578      END DO 
     1579 
     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 ) 
     1582 
     1583   END SUBROUTINE s_sh94 
     1584 
     1585   SUBROUTINE s_sf12 
     1586 
     1587      !!---------------------------------------------------------------------- 
     1588      !!                  ***  ROUTINE s_sf12 ***  
     1589      !!                      
     1590      !! ** Purpose :   stretch the s-coordinate system 
     1591      !! 
     1592      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012? 
     1593      !!                mixed S/sigma/Z coordinate 
     1594      !! 
     1595      !!                This method allows the maintenance of fixed surface and or 
     1596      !!                bottom cell resolutions (cf. geopotential coordinates)  
     1597      !!                within an analytically derived stretched S-coordinate framework. 
     1598      !! 
     1599      !! 
     1600      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 
     1601      !!---------------------------------------------------------------------- 
     1602      ! 
     1603      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     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 
     1618 
     1619      DO ji = 1, jpi 
     1620         DO jj = 1, jpj 
     1621 
     1622          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 
     1623               
     1624              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,. 
     1625                                                     ! could be changed by users but care must be taken to do so carefully 
     1626              zzb = 1.0_wp-(zzb/hbatt(ji,jj)) 
     1627             
     1628              zzs = rn_zs / hbatt(ji,jj)  
     1629               
     1630              IF (rn_efold /= 0.0_wp) THEN 
     1631                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 
     1632              ELSE 
     1633                zsmth = 1.0_wp  
     1634              ENDIF 
     1635                
     1636              DO jk = 1, jpk 
     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) 
     1639              ENDDO 
     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  ) 
     1642  
     1643          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma 
     1644 
     1645            DO jk = 1, jpk 
     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) 
     1648            END DO 
     1649 
     1650          ELSE  ! shallow water, z coordinates 
     1651 
     1652            DO jk = 1, jpk 
     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)) 
     1655            END DO 
     1656 
     1657          ENDIF 
     1658 
     1659          DO jk = 1, jpkm1 
     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) 
     1662          END DO 
     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)) 
     1665 
     1666          ! Coefficients for vertical depth as the sum of e3w scale factors 
     1667          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 
     1668          DO jk = 2, jpk 
     1669             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 
     1670          END DO 
     1671 
     1672          DO jk = 1, jpk 
     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) 
     1676          END DO 
     1677 
     1678        ENDDO   ! for all jj's 
     1679      ENDDO    ! for all ji's 
     1680 
     1681      DO ji=1,jpi 
     1682        DO jj=1,jpj 
     1683 
     1684          DO jk = 1, jpk 
     1685                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 
     1686                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1687                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 
     1688                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     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) ) / & 
     1691                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     1692                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 
     1693                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1694                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 
     1695                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1696 
     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) 
     1701             ! 
     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) 
     1705          END DO 
     1706 
     1707        ENDDO 
     1708      ENDDO 
     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 ) 
     1713 
     1714   END SUBROUTINE s_sf12 
     1715 
     1716   SUBROUTINE s_tanh() 
     1717 
     1718      !!---------------------------------------------------------------------- 
     1719      !!                  ***  ROUTINE s_tanh***  
     1720      !!                      
     1721      !! ** Purpose :   stretch the s-coordinate system 
     1722      !! 
     1723      !! ** Method  :   s-coordinate stretch  
     1724      !! 
     1725      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 
     1726      !!---------------------------------------------------------------------- 
     1727 
     1728      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     1729      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     1730 
     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 
    15771740      DO jk = 1, jpk 
    1578          DO jj = 1, jpj 
    1579             DO ji = 1, jpi 
    1580                IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
    1581                   WRITE(ctmp1,*) 'zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    1582                   CALL ctl_stop( ctmp1 ) 
    1583                ENDIF 
    1584                IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 
    1585                   WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
    1586                   CALL ctl_stop( ctmp1 ) 
    1587                ENDIF 
    1588             END DO 
    1589          END DO 
    1590       END DO 
    1591 !!gm bug    #endif 
    1592       ! 
    1593       CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           ) 
    1594       CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
    1595       CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
    1596       ! 
    1597       IF( nn_timing == 1 )  CALL timing_stop('zgr_sco') 
    1598       ! 
    1599    END SUBROUTINE zgr_sco 
     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) 
     1745      ! 
     1746      ! Coefficients for vertical scale factors at w-, t- levels 
     1747!!gm bug :  define it from analytical function, not like juste bellow.... 
     1748!!gm        or betteroffer the 2 possibilities.... 
     1749      DO jk = 1, jpkm1 
     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) ) 
     1755      ! 
     1756      ! Coefficients for vertical depth as the sum of e3w scale factors 
     1757      z_gsi3w(1) = 0.5_wp * z_esigw(1) 
     1758      DO jk = 2, jpk 
     1759         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk) 
     1760      END DO 
     1761!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
     1762      DO jk = 1, jpk 
     1763         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
     1764         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
     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 ) 
     1768      END DO 
     1769!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
     1770      DO jj = 1, jpj 
     1771         DO ji = 1, jpi 
     1772            DO jk = 1, jpk 
     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) ) 
     1777              ! 
     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 
     1788   END SUBROUTINE s_tanh 
     1789 
     1790   FUNCTION fssig( pk ) RESULT( pf ) 
     1791      !!---------------------------------------------------------------------- 
     1792      !!                 ***  ROUTINE fssig *** 
     1793      !!        
     1794      !! ** Purpose :   provide the analytical function in s-coordinate 
     1795      !!           
     1796      !! ** Method  :   the function provide the non-dimensional position of 
     1797      !!                T and W (i.e. between 0 and 1) 
     1798      !!                T-points at integer values (between 1 and jpk) 
     1799      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
     1800      !!---------------------------------------------------------------------- 
     1801      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate 
     1802      REAL(wp)             ::   pf   ! sigma value 
     1803      !!---------------------------------------------------------------------- 
     1804      ! 
     1805      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   & 
     1806         &     - TANH( rn_thetb * rn_theta                                )  )   & 
     1807         & * (   COSH( rn_theta                           )                      & 
     1808         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              & 
     1809         & / ( 2._wp * SINH( rn_theta ) ) 
     1810      ! 
     1811   END FUNCTION fssig 
     1812 
     1813 
     1814   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 
     1815      !!---------------------------------------------------------------------- 
     1816      !!                 ***  ROUTINE fssig1 *** 
     1817      !! 
     1818      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate 
     1819      !! 
     1820      !! ** Method  :   the function provides the non-dimensional position of 
     1821      !!                T and W (i.e. between 0 and 1) 
     1822      !!                T-points at integer values (between 1 and jpk) 
     1823      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
     1824      !!---------------------------------------------------------------------- 
     1825      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate 
     1826      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient 
     1827      REAL(wp)             ::   pf1   ! sigma value 
     1828      !!---------------------------------------------------------------------- 
     1829      ! 
     1830      IF ( rn_theta == 0 ) then      ! uniform sigma 
     1831         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 
     1832      ELSE                        ! stretched sigma 
     1833         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              & 
     1834            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  & 
     1835            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  ) 
     1836      ENDIF 
     1837      ! 
     1838   END FUNCTION fssig1 
     1839 
     1840 
     1841   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma ) 
     1842      !!---------------------------------------------------------------------- 
     1843      !!                 ***  ROUTINE fgamma  *** 
     1844      !! 
     1845      !! ** Purpose :   provide analytical function for the s-coordinate 
     1846      !! 
     1847      !! ** Method  :   the function provides the non-dimensional position of 
     1848      !!                T and W (i.e. between 0 and 1) 
     1849      !!                T-points at integer values (between 1 and jpk) 
     1850      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
     1851      !! 
     1852      !!                This method allows the maintenance of fixed surface and or 
     1853      !!                bottom cell resolutions (cf. geopotential coordinates)  
     1854      !!                within an analytically derived stretched S-coordinate framework. 
     1855      !! 
     1856      !! Reference  :   Siddorn and Furner, in prep 
     1857      !!---------------------------------------------------------------------- 
     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 
     1866      integer                 ::   jk 
     1867      !!---------------------------------------------------------------------- 
     1868      ! 
     1869 
     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) 
     1876      
     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 
     1881  
     1882      DO jk = 1, jpk 
     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) 
     1885      ENDDO  
     1886 
     1887      ! 
     1888   END FUNCTION fgamma 
    16001889 
    16011890   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.