Changeset 3453


Ignore:
Timestamp:
2012-08-14T15:38:51+02:00 (8 years ago)
Author:
rfurner
Message:

changes to allow new strecthed coordinate, including namelist for AMM12 std config; ticket#985

Location:
branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/AMM12/EXP00/namelist

    r3309 r3453  
    6565&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate 
    6666!----------------------------------------------------------------------- 
    67    rn_sbot_min =   10.     !  minimum depth of s-bottom surface (>0) (m) 
    68    rn_sbot_max = 7000.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     67   ln_s_sh94   = .false.   !  Song & Haidvogel 1994 hybrid S-sigma   (T)| 
     68   ln_s_sf12   = .true.    !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied 
     69   ln_sigcrit  = .true.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch  
     70   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m) 
     71   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     72   rn_hc       =   50.0    !  critical depth for transition to stretched coordinates 
     73   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1) 
     74                           !  SH94 stretching coefficients  
    6975   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20) 
    70    rn_thetb    =    1.00   !  bottom control parameter  (0<=thetb<= 1) 
    71    rn_rmax     =    0.30   !  maximum cut-off r-value allowed (0<r_max<1) 
    72    ln_s_sigma  = .true.    !  hybrid s-sigma coordinates 
    73    rn_bb       =    0.8    !  stretching with s-sigma 
    74    rn_hc       =  150.0    !  critical depth with s-sigma  
     76   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1) 
     77   rn_bb       =    0.8    !  stretching with SH94 s-sigma 
     78                           !  SF12 stretching coefficient  
     79   rn_alpha    =    4.4    !  stretching with SH94 s-sigma 
     80   rn_efold    =    0.0    !  efold length scale for transition to stretched coord 
     81   rn_zs       =    1.0    !  depth of surface grid box 
     82                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 
     83   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb 
     84   rn_zb_b     =   -0.2    !  offset for calculating Zb 
    7585/ 
    7686!----------------------------------------------------------------------- 
     
    7989   nn_bathy    =    1      !  compute (=0) or read (=1) the bathymetry file 
    8090   nn_closea    =   0      !  remove (=0) or keep (=1) closed seas and lakes (ORCA) 
    81    nn_msh      =    0      !  create (=1) a mesh file or not (=0) 
     91   nn_msh      =    1      !  create (=1) a mesh file or not (=0) 
    8292   rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0) 
    8393   rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of 
     
    419429   bn_v3d  =    'amm12_bdyV_u3d' ,         24        , 'vomecrty' ,     .true.     , .false. ,  'daily'  ,    ''    ,   '' 
    420430   bn_tem  =    'amm12_bdyT_tra' ,         24        , 'votemper' ,     .true.     , .false. ,  'daily'  ,    ''    ,   '' 
    421    bn_sal  =    'amm12_bdyT_tra' ,         24        , 'vosaline' ,     .true.     , .false. ,  'daily'  ,    ''    ,   '' 
     431   bn_tem  =    'amm12_bdyT_tra' ,         24        , 'vosaline' ,     .true.     , .false. ,  'daily'  ,    ''    ,   '' 
    422432   cn_dir  =    'bdydta/' 
    423433   ln_full_vel = .false. 
  • branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r3294 r3453  
    188188         CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 
    189189         ! 
    190          CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 )    !    ! stretched system 
    191          CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 ) 
     190         CALL iom_rstput( 0, 0, inum4, 'gdept' , gdept )    !    ! stretched system 
     191         CALL iom_rstput( 0, 0, inum4, 'gdepw' , gdepw )  
    192192      ENDIF 
    193193       
  • branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r3421 r3453  
    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 
    1718   !!---------------------------------------------------------------------- 
    1819 
     
    2829   !!       zgr_sco      : s-coordinate 
    2930   !!       fssig        : sigma coordinate non-dimensional function 
     31   !!       fgamma       : Siddorn and Furner stretching function 
    3032   !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!) 
    3133   !!--------------------------------------------------------------------- 
     
    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   LOGICAL  ::   ln_sigcrit  = .false.      ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch  
     54   ! 
    4955   REAL(wp) ::   rn_sbot_min =  300._wp     ! minimum depth of s-bottom surface (>0) (m) 
    5056   REAL(wp) ::   rn_sbot_max = 5250._wp     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     57   REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1) 
     58   REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for transition from sigma to stretched coordinates 
     59   ! Song and Haidvogel 1994 stretching parameters 
    5160   REAL(wp) ::   rn_theta    =    6.00_wp   ! surface control parameter (0<=rn_theta<=20) 
    5261   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 
     62   REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter  
    5663   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 
    57    REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for s-sigma coordinates 
     64   ! Siddorn and Furner stretching parameters 
     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      !!---------------------------------------------------------------------- 
     
    11101071      !!            esigt(k) = fsdsig(k    ) 
    11111072      !!            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 
     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,*) '        Siddorn and Furner 2012 stretching coefficients' 
     1128         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha 
     1129         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold 
     1130         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs 
     1131         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a 
     1132         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b 
     1133         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 
     1134      ENDIF 
    11611135 
    11621136      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate 
     
    13521326      ! non-dimensional "sigma" for model level depth at w- and t-levels 
    13531327 
    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  
     1328 
     1329!======================================================================== 
     1330! Song and Haidvogel  1994 (ln_s_sh94=T) 
     1331! Siddorn and Furner 2012 (ln_sf12=T) 
     1332! or  tanh function       (both false)                     
     1333!======================================================================== 
     1334      IF      ( ln_s_sh94 ) THEN  
     1335                           CALL s_sh94() 
     1336      ELSE IF ( ln_s_sf12 ) THEN 
     1337                           CALL s_sf12() 
     1338      ELSE                  
     1339                           CALL s_tanh() 
     1340      ENDIF  
     1341 
     1342      CALL lbc_lnk( e3t , 'T', 1._wp ) 
     1343      CALL lbc_lnk( e3u , 'U', 1._wp ) 
     1344      CALL lbc_lnk( e3v , 'V', 1._wp ) 
     1345      CALL lbc_lnk( e3f , 'F', 1._wp ) 
     1346      CALL lbc_lnk( e3w , 'W', 1._wp ) 
     1347      CALL lbc_lnk( e3uw, 'U', 1._wp ) 
     1348      CALL lbc_lnk( e3vw, 'V', 1._wp ) 
     1349 
     1350      fsdepw(:,:,:) = gdepw (:,:,:) 
     1351      fsde3w(:,:,:) = gdep3w(:,:,:) 
    14891352      ! 
    14901353      where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1.0 
     
    15441407            &                          ' w ', MAXVAL( fse3w (:,:,:) ) 
    15451408      ENDIF 
    1546       ! 
     1409      !  END DO 
    15471410      IF(lwp) THEN                                  ! selected vertical profiles 
    15481411         WRITE(numout,*) 
     
    15741437      ENDIF 
    15751438 
    1576 !!gm bug?  no more necessary?  if ! defined key_helsinki 
    1577       DO jk = 1, jpk 
     1439!================================================================================ 
     1440! check the coordinate makes sense 
     1441!================================================================================ 
     1442      DO ji = 1, jpi 
    15781443         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 
     1444 
     1445            IF( hbatt(ji,jj) > 0._wp) THEN 
     1446               DO jk = 1, mbathy(ji,jj) 
     1447                 ! check coordinate is monotonically increasing 
     1448                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
     1449                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1450                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1451                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 
     1452                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 
     1453                    CALL ctl_stop( ctmp1 ) 
     1454                 ENDIF 
     1455                 ! and check it has never gone negative 
     1456                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 
     1457                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
     1458                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1459                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
     1460                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
     1461                    CALL ctl_stop( ctmp1 ) 
     1462                 ENDIF 
     1463                 ! and check it never exceeds the total depth 
     1464                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1465                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1466                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1467                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
     1468                    CALL ctl_stop( ctmp1 ) 
     1469                 ENDIF 
     1470               END DO 
     1471 
     1472               DO jk = 1, mbathy(ji,jj)-1 
     1473                 ! and check it never exceeds the total depth 
     1474                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1475                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1476                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1477                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
     1478                    CALL ctl_stop( ctmp1 ) 
     1479                 ENDIF 
     1480               END DO 
     1481 
     1482            ENDIF 
     1483 
     1484         END DO 
     1485      END DO 
    15921486      ! 
    15931487      CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           ) 
     1488      ! 
     1489      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco') 
     1490      ! 
     1491   END SUBROUTINE zgr_sco 
     1492 
     1493!!====================================================================== 
     1494   SUBROUTINE s_sh94() 
     1495 
     1496      !!---------------------------------------------------------------------- 
     1497      !!                  ***  ROUTINE s_sh94  *** 
     1498      !!                      
     1499      !! ** Purpose :   stretch the s-coordinate system 
     1500      !! 
     1501      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994 
     1502      !!                mixed S/sigma coordinate 
     1503      !! 
     1504      !! Reference : Song and Haidvogel 1994.  
     1505      !!---------------------------------------------------------------------- 
     1506      ! 
     1507      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     1508      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     1509      ! 
     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 
     1520 
     1521      DO ji = 1, jpi 
     1522         DO jj = 1, jpj 
     1523 
     1524            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
     1525               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 ) 
     1528               END DO 
     1529            ELSE ! shallow water, uniform sigma 
     1530               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) 
     1533                  END DO 
     1534            ENDIF 
     1535            ! 
     1536            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) ) 
     1542            ! 
     1543            ! Coefficients for vertical depth as the sum of e3w scale factors 
     1544            gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 
     1545            DO jk = 2, jpk 
     1546               gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
     1547            END DO 
     1548            ! 
     1549            DO jk = 1, jpk 
     1550               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
     1551               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 ) 
     1555            END DO 
     1556           ! 
     1557         END DO   ! for all jj's 
     1558      END DO    ! for all ji's 
     1559 
     1560      DO ji = 1, jpim1 
     1561         DO jj = 1, jpjm1 
     1562            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) )   & 
     1564                  &              / ( 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) )   & 
     1566                  &              / ( 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) )   & 
     1569                  &              / ( 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) )   & 
     1571                  &              / ( 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) )   & 
     1573                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1574               ! 
     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) ) 
     1579               ! 
     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) ) 
     1583            END DO 
     1584        END DO 
     1585      END DO 
     1586 
    15941587      CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
    15951588      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 
     1589 
     1590   END SUBROUTINE s_sh94 
     1591 
     1592   SUBROUTINE s_sf12 
     1593 
     1594      !!---------------------------------------------------------------------- 
     1595      !!                  ***  ROUTINE s_sf12 ***  
     1596      !!                      
     1597      !! ** Purpose :   stretch the s-coordinate system 
     1598      !! 
     1599      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012? 
     1600      !!                mixed S/sigma/Z coordinate 
     1601      !! 
     1602      !!                This method allows the maintenance of fixed surface and or 
     1603      !!                bottom cell resolutions (cf. geopotential coordinates)  
     1604      !!                within an analytically derived stretched S-coordinate framework. 
     1605      !! 
     1606      !! 
     1607      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 
     1608      !!---------------------------------------------------------------------- 
     1609      ! 
     1610      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 
     1625 
     1626      DO ji = 1, jpi 
     1627         DO jj = 1, jpj 
     1628 
     1629          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 
     1630               
     1631              zbb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,. 
     1632                                                     ! could be changed by users but care must be taken to do so carefully 
     1633              zbb = 1.0_wp-(zbb/hbatt(ji,jj)) 
     1634             
     1635              zss = rn_zs / hbatt(ji,jj)  
     1636               
     1637              IF (rn_efold /= 0.0_wp) THEN 
     1638                fsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 
     1639              ELSE 
     1640                fsmth = 1.0_wp  
     1641              ENDIF 
     1642                
     1643              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) 
     1646              ENDDO 
     1647              gsigw3(ji,jj,:) = fgamma( gsigw3(ji,jj,:), zbb, zss, fsmth  ) 
     1648              gsigt3(ji,jj,:) = fgamma( gsigt3(ji,jj,:), zbb, zss, fsmth  ) 
     1649  
     1650          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma 
     1651 
     1652            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) 
     1655            END DO 
     1656 
     1657          ELSE  ! shallow water, z coordinates 
     1658 
     1659            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)) 
     1662            END DO 
     1663 
     1664          ENDIF 
     1665 
     1666          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) 
     1669          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)) 
     1672 
     1673          ! Coefficients for vertical depth as the sum of e3w scale factors 
     1674          gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1) 
     1675          DO jk = 2, jpk 
     1676             gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
     1677          END DO 
     1678 
     1679          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) 
     1683          END DO 
     1684 
     1685        ENDDO   ! for all jj's 
     1686      ENDDO    ! for all ji's 
     1687 
     1688      DO ji=1,jpi 
     1689        DO jj=1,jpj 
     1690 
     1691          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) ) / & 
     1693                                    ( 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) ) / & 
     1695                                    ( 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) ) / & 
     1698                                    ( 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) ) / & 
     1700                                    ( 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) ) / & 
     1702                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1703 
     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) 
     1708             ! 
     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) 
     1712          END DO 
     1713 
     1714        ENDDO 
     1715      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 ) 
     1719 
     1720   END SUBROUTINE s_sf12 
     1721 
     1722   SUBROUTINE s_tanh() 
     1723 
     1724      !!---------------------------------------------------------------------- 
     1725      !!                  ***  ROUTINE s_tanh***  
     1726      !!                      
     1727      !! ** Purpose :   stretch the s-coordinate system 
     1728      !! 
     1729      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012? 
     1730      !!                mixed S/sigma/Z coordinate 
     1731      !! 
     1732      !! 
     1733      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 
     1734      !!---------------------------------------------------------------------- 
     1735 
     1736      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     1737      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     1738 
     1739      DO jk = 1, jpk 
     1740        gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 
     1741        gsigt(jk) = -fssig( REAL(jk,wp)        ) 
     1742      END DO 
     1743      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk) 
     1744      ! 
     1745      ! Coefficients for vertical scale factors at w-, t- levels 
     1746!!gm bug :  define it from analytical function, not like juste bellow.... 
     1747!!gm        or betteroffer the 2 possibilities.... 
     1748      DO jk = 1, jpkm1 
     1749         esigt(jk  ) = gsigw(jk+1) - gsigw(jk) 
     1750         esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 
     1751      END DO 
     1752      esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) )  
     1753      esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 
     1754      ! 
     1755      ! Coefficients for vertical depth as the sum of e3w scale factors 
     1756      gsi3w(1) = 0.5_wp * esigw(1) 
     1757      DO jk = 2, jpk 
     1758         gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 
     1759      END DO 
     1760!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
     1761      DO jk = 1, jpk 
     1762         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
     1763         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
     1764         gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 
     1765         gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 
     1766         gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 
     1767      END DO 
     1768!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
     1769      DO jj = 1, jpj 
     1770         DO ji = 1, jpi 
     1771            DO jk = 1, jpk 
     1772              e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
     1773              e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
     1774              e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
     1775              e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 
     1776              ! 
     1777              e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
     1778              e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
     1779              e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
     1780            END DO 
     1781         END DO 
     1782      END DO 
     1783   END SUBROUTINE s_tanh 
     1784 
     1785   FUNCTION fssig( pk ) RESULT( pf ) 
     1786      !!---------------------------------------------------------------------- 
     1787      !!                 ***  ROUTINE fssig *** 
     1788      !!        
     1789      !! ** Purpose :   provide the analytical function in s-coordinate 
     1790      !!           
     1791      !! ** Method  :   the function provide the non-dimensional position of 
     1792      !!                T and W (i.e. between 0 and 1) 
     1793      !!                T-points at integer values (between 1 and jpk) 
     1794      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
     1795      !!---------------------------------------------------------------------- 
     1796      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate 
     1797      REAL(wp)             ::   pf   ! sigma value 
     1798      !!---------------------------------------------------------------------- 
     1799      ! 
     1800      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   & 
     1801         &     - TANH( rn_thetb * rn_theta                                )  )   & 
     1802         & * (   COSH( rn_theta                           )                      & 
     1803         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              & 
     1804         & / ( 2._wp * SINH( rn_theta ) ) 
     1805      ! 
     1806   END FUNCTION fssig 
     1807 
     1808 
     1809   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 
     1810      !!---------------------------------------------------------------------- 
     1811      !!                 ***  ROUTINE fssig1 *** 
     1812      !! 
     1813      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate 
     1814      !! 
     1815      !! ** Method  :   the function provides the non-dimensional position of 
     1816      !!                T and W (i.e. between 0 and 1) 
     1817      !!                T-points at integer values (between 1 and jpk) 
     1818      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
     1819      !!---------------------------------------------------------------------- 
     1820      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate 
     1821      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient 
     1822      REAL(wp)             ::   pf1   ! sigma value 
     1823      !!---------------------------------------------------------------------- 
     1824      ! 
     1825      IF ( rn_theta == 0 ) then      ! uniform sigma 
     1826         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 
     1827      ELSE                        ! stretched sigma 
     1828         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              & 
     1829            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  & 
     1830            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  ) 
     1831      ENDIF 
     1832      ! 
     1833   END FUNCTION fssig1 
     1834 
     1835 
     1836   FUNCTION fgamma( pk1, Zbb, Zss, F ) RESULT( gam ) 
     1837      !!---------------------------------------------------------------------- 
     1838      !!                 ***  ROUTINE fgamma  *** 
     1839      !! 
     1840      !! ** Purpose :   provide analytical function for the s-coordinate 
     1841      !! 
     1842      !! ** Method  :   the function provides the non-dimensional position of 
     1843      !!                T and W (i.e. between 0 and 1) 
     1844      !!                T-points at integer values (between 1 and jpk) 
     1845      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
     1846      !! 
     1847      !!                This method allows the maintenance of fixed surface and or 
     1848      !!                bottom cell resolutions (cf. geopotential coordinates)  
     1849      !!                within an analytically derived stretched S-coordinate framework. 
     1850      !! 
     1851      !! Reference  :   Siddorn and Furner, in prep 
     1852      !!---------------------------------------------------------------------- 
     1853      REAL(wp), INTENT(in   ) ::   pk1(jpk)   ! continuous "k" coordinate 
     1854      REAL(wp)                ::   gam(jpk)   ! stretched coordinate 
     1855      REAL(wp), INTENT(in   ) ::   Zbb      ! Bottom box depth 
     1856      REAL(wp), INTENT(in   ) ::   Zss      ! surface box depth 
     1857      REAL(wp), INTENT(in   ) ::   F        ! Smoothing parameter 
     1858      REAL(wp)                ::   a1,a2,a3 ! local variables 
     1859      REAL(wp)                ::   n1,n2    ! local variables 
     1860      REAL(wp)                ::   A,B,X    ! local variables 
     1861      integer                 ::   jk 
     1862      !!---------------------------------------------------------------------- 
     1863      ! 
     1864 
     1865      n1  =  1./(jpk-1.) 
     1866      n2  =  1. -  n1 
     1867 
     1868      a1 = (rn_alpha+2.0_wp)*n1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n1**(rn_alpha+2.0_wp)  
     1869      a2 = (rn_alpha+2.0_wp)*n2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n2**(rn_alpha+2.0_wp) 
     1870      a3 = ( n2**3.0_wp - a2)/( n1**3.0_wp - a1) 
     1871      
     1872      A = Zbb - a3*(Zss-a1)-a2 
     1873      A = A/( n2-0.5_wp*(a2+n2**2.0_wp) - a3*(n1-0.5_wp*(a1+n1**2.0_wp) ) ) 
     1874      B = (Zss - a1 - A*( n1-0.5_wp*(a1+n1**2.0_wp ) ) ) / (n1**3.0_wp - a1) 
     1875      X = 1.0_wp-A/2.0_wp-B 
     1876  
     1877      DO jk = 1, jpk 
     1878        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) ) 
     1879        gam(jk) = gam(jk)*F+pk1(jk)*(1.0_wp-F) 
     1880      ENDDO  
     1881 
     1882      ! 
     1883   END FUNCTION fgamma 
    16001884 
    16011885   !!====================================================================== 
  • branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/par_AMM_12km.h90

    r3294 r3453  
    1919      jpidta  = 198,        &  !: first horizontal dimension > or = to jpi 
    2020      jpjdta  = 224,        &  !: second                     > or = to jpj 
    21       jpkdta  = 33,         &  !: number of levels           > or = to jpk 
     21      jpkdta  = 51,         &  !: number of levels           > or = to jpk 
    2222      ! total domain matrix size 
    2323      jpiglo  = jpidta,      &  !: first  dimension of global domain --> i 
Note: See TracChangeset for help on using the changeset viewer.