Changeset 5257


Ignore:
Timestamp:
2015-05-08T16:18:50+02:00 (5 years ago)
Author:
timgraham
Message:

Fixed lots of compilation errors

Location:
branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90

    r5255 r5257  
    4444      !!---------------------------------------------------------------------- 
    4545      ! 
     46      USE utils 
    4647      IMPLICIT NONE 
    4748 
    48       USE utils 
    4949 
    5050     !!---------------------------------------------------------------------- 
     
    5454      READ( numnam, namzgr_sco ) 
    5555      CLOSE( numnam ) 
    56       jpk = rn_jpk 
    57  
    58       WRITE(numout,*) 
    59       WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate' 
    60       WRITE(numout,*) '~~~~~~~~~~~' 
    61       WRITE(numout,*) '   Namelist namzgr_sco' 
    62       WRITE(numout,*) '     stretching coeffs ' 
    63       WRITE(numout,*) '        Number of levels                             rn_jpk        = ',rn_jpk 
    64       WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max 
    65       WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min 
    66       WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc 
    67       WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax 
    68       WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94 
    69       WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients' 
    70       WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta 
    71       WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb 
    72       WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb 
    73       WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12 
    74       WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit 
    75       WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients' 
    76       WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha 
    77       WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold 
    78       WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs 
    79       WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a 
    80       WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b 
    81       WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 
     56      jpk = INT(rn_jpk) 
     57 
     58      WRITE(*,*) 
     59      WRITE(*,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate' 
     60      WRITE(*,*) '~~~~~~~~~~~' 
     61      WRITE(*,*) '   Namelist namzgr_sco' 
     62      WRITE(*,*) '     stretching coeffs ' 
     63      WRITE(*,*) '        Number of levels                             rn_jpk        = ',rn_jpk 
     64      WRITE(*,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max 
     65      WRITE(*,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min 
     66      WRITE(*,*) '        Critical depth                               rn_hc         = ',rn_hc 
     67      WRITE(*,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax 
     68      WRITE(*,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94 
     69      WRITE(*,*) '        Song and Haidvogel 1994 stretching coefficients' 
     70      WRITE(*,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta 
     71      WRITE(*,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb 
     72      WRITE(*,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb 
     73      WRITE(*,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12 
     74      WRITE(*,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit 
     75      WRITE(*,*) '        Siddorn and Furner 2012 stretching coefficients' 
     76      WRITE(*,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha 
     77      WRITE(*,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold 
     78      WRITE(*,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs 
     79      WRITE(*,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a 
     80      WRITE(*,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b 
     81      WRITE(*,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 
    8282 
    8383       
     
    102102      DO jj = 1, jpj 
    103103         DO ji = 1, jpi 
    104            IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 
     104           IF( bathy(ji,jj) > 0. )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 
    105105         END DO 
    106106      END DO 
     
    114114      DO jj = 1, jpj 
    115115         DO ji = 1, jpi 
    116            IF( bathy(ji,jj) == 0._wp ) THEN 
     116           IF( bathy(ji,jj) == 0. ) THEN 
    117117             iip1 = MIN( ji+1, jpi ) 
    118118             ijp1 = MIN( jj+1, jpj ) 
     
    120120             ijm1 = MAX( jj-1, 1 ) 
    121121             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              & 
    122         &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 
     122        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0. ) THEN 
    123123               zenv(ji,jj) = rn_sbot_min 
    124124             ENDIF 
     
    128128      !  
    129129      ! smooth the bathymetry (if required) 
    130       scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea) 
     130      scosrf(:,:) = 0.             ! ocean surface depth (here zero: no under ice-shelf sea) 
    131131      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth 
    132132      ! 
    133133      jl = 0 
    134       zrmax = 1._wp 
     134      zrmax = 1. 
    135135      !    
    136136      !      
    137137      ! set scaling factor used in reducing vertical gradients 
    138       zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax ) 
     138      zrfact = ( 1. - rn_rmax ) / ( 1. + rn_rmax ) 
    139139      ! 
    140140      ! initialise temporary evelope depth arrays 
     
    145145      ! 
    146146      ! initialise temporary r-value arrays 
    147       zri(:,:) = 1._wp 
    148       zrj(:,:) = 1._wp 
     147      zri(:,:) = 1. 
     148      zrj(:,:) = 1. 
    149149      !                                                            ! ================ ! 
    150       DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  ! 
     150      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8 ) !  Iterative loop  ! 
    151151         !                                                         ! ================ ! 
    152152         jl = jl + 1 
    153          zrmax = 0._wp 
     153         zrmax = 0. 
    154154         ! we set zrmax from previous r-values (zri and zrj) first 
    155155         ! if set after current r-value calculation (as previously) 
    156156         ! we could exit DO WHILE prematurely before checking r-value 
    157157         ! of current zenv 
    158          DO jj = 1, nlcj 
    159             DO ji = 1, nlci 
     158!         DO jj = 1, nlcj 
     159!            DO ji = 1, nlci 
     160         DO jj = 1, jpi !jpi or jpim1? 
     161            DO ji = 1, jpj 
    160162               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) 
    161163            END DO 
    162164         END DO 
    163          zri(:,:) = 0._wp 
    164          zrj(:,:) = 0._wp 
    165          DO jj = 1, nlcj 
    166             DO ji = 1, nlci 
    167                iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi) 
    168                ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj) 
    169                IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 
     165         zri(:,:) = 0. 
     166         zrj(:,:) = 0. 
     167!         DO jj = 1, nlci 
     168 !           DO ji = 1, nlcj 
     169 !              iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi) 
     170 !              ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj) 
     171 !              IF( (zenv(ji,jj) > 0.) .AND. (zenv(iip1,jj) > 0.)) THEN 
     172 !                 zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) ) 
     173 !              END IF 
     174 !              IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) THEN 
     175 !                 zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) ) 
     176 !              END IF 
     177 !              IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact 
     178 !              IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact 
     179 !              IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact 
     180 !              IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact 
     181 !           END DO 
     182 !        END DO 
     183         DO jj = 1, jpi-1 
     184            DO ji = 1, jpj-1 
     185               iip1 = ji+1       
     186               ijp1 = jj+1       
     187               IF( (zenv(ji,jj) > 0.) .AND. (zenv(iip1,jj) > 0.)) THEN 
    170188                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) ) 
    171189               END IF 
    172                IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 
     190               IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) THEN 
    173191                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) ) 
    174192               END IF 
     
    182200         WRITE(*,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax 
    183201         ! 
    184          DO jj = 1, nlcj 
    185             DO ji = 1, nlci 
     202         DO jj = 1, jpi 
     203            DO ji = 1, jpj 
    186204               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) 
    187205            END DO 
    188206         END DO 
    189          ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero 
    190207         !                                                  ! ================ ! 
    191208      END DO                                                !     End loop     ! 
     
    198215      ! 
    199216      ! Envelope bathymetry saved in hbatt 
     217      ! TODO - get this section to work 
    200218      hbatt(:,:) = zenv(:,:)  
    201       IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 
    202          CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 
    203          DO jj = 1, jpj 
    204             DO ji = 1, jpi 
    205                ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp ) 
    206                hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper ) 
    207             END DO 
    208          END DO 
    209       ENDIF 
     219!      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0. ) THEN 
     220 !        CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 
     221 !        DO jj = 1, jpj 
     222 !           DO ji = 1, jpi 
     223 !              ztaper = EXP( -(gphit(ji,jj)/8.)**2. ) 
     224 !              hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1. - ztaper ) 
     225 !           END DO 
     226 !        END DO 
     227 !     ENDIF 
    210228      ! 
    211229      !                                        ! ============================== 
    212230      !                                        !   hbatu, hbatv, hbatf fields 
    213231      !                                        ! ============================== 
    214       WRITE(numout,*) 
    215       WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
     232      WRITE(*,*) 
     233      WRITE(*,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
    216234      hbatu(:,:) = rn_sbot_min 
    217235      hbatv(:,:) = rn_sbot_min 
    218236      hbatf(:,:) = rn_sbot_min 
    219       DO jj = 1, jpjm1 
    220         DO ji = 1, jpim1   ! NO vector opt. 
    221            hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) ) 
    222            hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) ) 
    223            hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   & 
     237      DO jj = 1, jpj-1 
     238        DO ji = 1, jpi-1   ! NO vector opt. 
     239           hbatu(ji,jj) = 0.50 * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) ) 
     240           hbatv(ji,jj) = 0.50 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) ) 
     241           hbatf(ji,jj) = 0.25 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   & 
    224242              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 
    225243        END DO 
     
    229247      DO jj = 1, jpj 
    230248         DO ji = 1, jpi 
    231             IF( hbatu(ji,jj) == 0._wp ) THEN 
    232                IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min 
    233                IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj) 
     249            IF( hbatu(ji,jj) == 0. ) THEN 
     250               IF( zhbat(ji,jj) == 0. )   hbatu(ji,jj) = rn_sbot_min 
     251               IF( zhbat(ji,jj) /= 0. )   hbatu(ji,jj) = zhbat(ji,jj) 
    234252            ENDIF 
    235253         END DO 
     
    238256      DO jj = 1, jpj 
    239257         DO ji = 1, jpi 
    240             IF( hbatv(ji,jj) == 0._wp ) THEN 
    241                IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min 
    242                IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj) 
     258            IF( hbatv(ji,jj) == 0. ) THEN 
     259               IF( zhbat(ji,jj) == 0. )   hbatv(ji,jj) = rn_sbot_min 
     260               IF( zhbat(ji,jj) /= 0. )   hbatv(ji,jj) = zhbat(ji,jj) 
    243261            ENDIF 
    244262         END DO 
     
    247265      DO jj = 1, jpj 
    248266         DO ji = 1, jpi 
    249             IF( hbatf(ji,jj) == 0._wp ) THEN 
    250                IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min 
    251                IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj) 
     267            IF( hbatf(ji,jj) == 0. ) THEN 
     268               IF( zhbat(ji,jj) == 0. )   hbatf(ji,jj) = rn_sbot_min 
     269               IF( zhbat(ji,jj) /= 0. )   hbatf(ji,jj) = zhbat(ji,jj) 
    252270            ENDIF 
    253271         END DO 
     
    260278      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
    261279 
    262          WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  & 
     280         WRITE(*,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  & 
    263281            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) 
    264          WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  & 
     282         WRITE(*,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  & 
    265283            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) ) 
    266          WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  & 
     284         WRITE(*,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  & 
    267285            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) 
    268          WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  & 
     286         WRITE(*,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  & 
    269287            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) 
    270288!! helsinki 
     
    326344         DO ji = 1, jpi 
    327345            DO jk = 1, jpk-1 
    328                IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    329             END DO 
    330             IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0 
    331          END DO 
    332       END DO 
    333       WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 
     346               IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
     347            END DO 
     348            IF( scobot(ji,jj) == 0.               )   mbathy(ji,jj) = 0 
     349         END DO 
     350      END DO 
     351      WRITE(*,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 
    334352 
    335353      ! min max values over the domain 
    336       WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) ) 
    337       WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   & 
     354      WRITE(*,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) ) 
     355      WRITE(*,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   & 
    338356         &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) ) 
    339       WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   & 
     357      WRITE(*,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   & 
    340358         &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   & 
    341359         &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   & 
    342360         &                          ' w ', MINVAL( e3w_0  (:,:,:) ) 
    343361 
    344       WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   & 
     362      WRITE(*,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   & 
    345363         &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) ) 
    346       WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   & 
     364      WRITE(*,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   & 
    347365         &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   & 
    348366         &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   & 
     
    350368       
    351369         ! selected vertical profiles 
    352       WRITE(numout,*) 
    353       WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 
    354       WRITE(numout,*) ' ~~~~~~  --------------------' 
    355       WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
    356       WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     & 
     370      WRITE(*,*) 
     371      WRITE(*,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 
     372      WRITE(*,*) ' ~~~~~~  --------------------' 
     373      WRITE(*,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
     374      WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     & 
    357375         &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk ) 
    358       DO jj = mj0(20), mj1(20) 
    359          DO ji = mi0(20), mi1(20) 
    360             WRITE(numout,*) 
    361             WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
    362             WRITE(numout,*) ' ~~~~~~  --------------------' 
    363             WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
    364             WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     & 
     376      DO jj = 20, 20 
     377         DO ji = 20, 20 
     378            WRITE(*,*) 
     379            WRITE(*,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
     380            WRITE(*,*) ' ~~~~~~  --------------------' 
     381            WRITE(*,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
     382            WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     & 
    365383               &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 
    366384         END DO 
    367385      END DO 
    368       DO jj = mj0(74), mj1(74) 
    369          DO ji = mi0(100), mi1(100) 
    370             WRITE(numout,*) 
    371             WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
    372             WRITE(numout,*) ' ~~~~~~  --------------------' 
    373             WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
    374             WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     & 
     386      DO jj = 74,74 
     387         DO ji = 100, 100 
     388            WRITE(*,*) 
     389            WRITE(*,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
     390            WRITE(*,*) ' ~~~~~~  --------------------' 
     391            WRITE(*,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
     392            WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     & 
    375393               &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 
    376394         END DO 
     
    383401         DO jj = 1, jpj 
    384402 
    385             IF( hbatt(ji,jj) > 0._wp) THEN 
     403            IF( hbatt(ji,jj) > 0.) THEN 
    386404               DO jk = 1, mbathy(ji,jj) 
    387405                 ! check coordinate is monotonically increasing 
    388                  IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
    389                     WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    390                     WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 
    391                     WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 
    392                     CALL STOP 1 
     406                 IF (e3w_0(ji,jj,jk) <= 0. .OR. e3t_0(ji,jj,jk) <= 0. ) THEN 
     407                    WRITE(*,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     408                    WRITE(*,*) 'e3w',e3w_0(ji,jj,:) 
     409                    WRITE(*,*) 'e3t',e3t_0(ji,jj,:) 
     410                    STOP 1 
    393411                 ENDIF 
    394412                 ! and check it has never gone negative 
    395                  IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 
    396                     WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
    397                     WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
    398                     WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
    399                     CALL STOP 1 
     413                 IF( gdepw_0(ji,jj,jk) < 0. .OR. gdept_0(ji,jj,jk) < 0. ) THEN 
     414                    WRITE(*,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
     415                    WRITE(*,*) 'gdepw',gdepw_0(ji,jj,:) 
     416                    WRITE(*,*) 'gdept',gdept_0(ji,jj,:) 
     417                    STOP 1 
    400418                 ENDIF 
    401419                 ! and check it never exceeds the total depth 
    402                  IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    403                     WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
    404                     WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
    405                     CALL STOP 1 
     420                 IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     421                    WRITE(*,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
     422                    WRITE(*,*) 'gdepw',gdepw_0(ji,jj,:) 
     423                    STOP 1 
    406424                 ENDIF 
    407425               END DO 
     
    409427               DO jk = 1, mbathy(ji,jj)-1 
    410428                 ! and check it never exceeds the total depth 
    411                 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    412                     WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
    413                     WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
    414                     CALL STOP 1 
     429                IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     430                    WRITE(*,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
     431                    WRITE(*,*) 'gdept',gdept_0(ji,jj,:) 
     432                    STOP 1 
    415433                 ENDIF 
    416434               END DO 
     
    421439      END DO 
    422440      ! 
    423       ! 
     441      ! Write output file 
     442      CALL write_coord_file() 
     443      ! 
     444      ! 
     445END PROGRAM SCOORD_GEN 
    424446 
    425447!!====================================================================== 
     
    437459      !!---------------------------------------------------------------------- 
    438460      ! 
    439       INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     461      USE utils 
    440462      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
    441463      ! 
     
    449471      ALLOCATE( z_esigwv3(jpi,jpj,jpk) ) 
    450472 
    451       z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp 
    452       z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp  
    453       z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp 
    454       z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp 
     473      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0. 
     474      z_esigt3  = 0.   ;   z_esigw3  = 0.  
     475      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0. 
     476      z_esigwu3 = 0.   ;   z_esigwv3 = 0. 
    455477 
    456478      DO ji = 1, jpi 
     
    459481            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
    460482               DO jk = 1, jpk 
    461                   z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 
     483                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5, rn_bb ) 
    462484                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb ) 
    463485               END DO 
     
    465487               DO jk = 1, jpk 
    466488                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp) 
    467                   z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 
     489                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp) 
    468490                  END DO 
    469491            ENDIF 
     
    473495               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 
    474496            END DO 
    475             z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) ) 
    476             z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) 
     497            z_esigw3(ji,jj,1  ) = 2. * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) ) 
     498            z_esigt3(ji,jj,jpk) = 2. * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) 
    477499            ! 
    478500            ! Coefficients for vertical depth as the sum of e3w scale factors 
    479             z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1) 
     501            z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 
    480502            DO jk = 2, jpk 
    481503               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 
     
    483505            ! 
    484506            DO jk = 1, jpk 
    485                zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpk-1,wp) 
    486                zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpk-1,wp) 
    487                gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
    488                gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
     507               zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) 
     508               zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) 
     509               gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
     510               gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
    489511               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 
    490512            END DO 
     
    543565      !!---------------------------------------------------------------------- 
    544566      ! 
    545       INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     567      USE utils 
    546568      REAL(wp) ::   zsmth               ! smoothing around critical depth 
    547569      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space 
     
    556578      ALLOCATE( z_esigwv3(jpi,jpj,jpk) ) 
    557579 
    558       z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp 
    559       z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp  
    560       z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp 
    561       z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp 
     580      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0. 
     581      z_esigt3  = 0.   ;   z_esigw3  = 0.  
     582      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0. 
     583      z_esigwu3 = 0.   ;   z_esigwv3 = 0. 
    562584 
    563585      DO ji = 1, jpi 
     
    568590              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,. 
    569591                                                     ! could be changed by users but care must be taken to do so carefully 
    570               zzb = 1.0_wp-(zzb/hbatt(ji,jj)) 
     592              zzb = 1.0-(zzb/hbatt(ji,jj)) 
    571593             
    572594              zzs = rn_zs / hbatt(ji,jj)  
    573595               
    574               IF (rn_efold /= 0.0_wp) THEN 
     596              IF (rn_efold /= 0.0) THEN 
    575597                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 
    576598              ELSE 
    577                 zsmth = 1.0_wp  
     599                zsmth = 1.0  
    578600              ENDIF 
    579601                
    580602              DO jk = 1, jpk 
    581603                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp) 
    582                 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp) 
     604                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 
    583605              ENDDO 
    584606              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  ) 
    585607              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  ) 
    586608  
    587           ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma 
     609          ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma 
    588610 
    589611            DO jk = 1, jpk 
     
    596618            DO jk = 1, jpk 
    597619              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
    598               z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
     620              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 
    599621            END DO 
    600622 
     
    605627             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 
    606628          END DO 
    607           z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  )) 
    608           z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) 
     629          z_esigw3(ji,jj,1  ) = 2.0 * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  )) 
     630          z_esigt3(ji,jj,jpk) = 2.0 * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) 
    609631 
    610632          ! Coefficients for vertical depth as the sum of e3w scale factors 
     
    615637 
    616638          DO jk = 1, jpk 
    617              gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 
    618              gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 
     639             gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 
     640             gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 
    619641             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 
    620642          END DO 
     
    671693      !!---------------------------------------------------------------------- 
    672694 
    673       INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     695      USE utils 
    674696      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
    675697 
     
    680702      ALLOCATE( z_esigt(jpk), z_esigw(jpk) ) 
    681703 
    682       z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp 
    683       z_esigt  = 0._wp   ;   z_esigw  = 0._wp  
     704      z_gsigw  = 0.   ;   z_gsigt  = 0.   ;   z_gsi3w  = 0. 
     705      z_esigt  = 0.   ;   z_esigw  = 0.  
    684706 
    685707      DO jk = 1, jpk 
    686         z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 
     708        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5 ) 
    687709        z_gsigt(jk) = -fssig( REAL(jk,wp)        ) 
    688710      END DO 
     
    696718         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk) 
    697719      END DO 
    698       z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) )  
    699       z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) ) 
     720      z_esigw( 1 ) = 2. * ( z_gsigt(1  ) - z_gsigw(1  ) )  
     721      z_esigt(jpk) = 2. * ( z_gsigt(jpk) - z_gsigw(jpk) ) 
    700722      ! 
    701723      ! Coefficients for vertical depth as the sum of e3w scale factors 
    702       z_gsi3w(1) = 0.5_wp * z_esigw(1) 
     724      z_gsi3w(1) = 0.5 * z_esigw(1) 
    703725      DO jk = 2, jpk 
    704726         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk) 
     
    706728!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
    707729      DO jk = 1, jpk 
    708          zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpk-1,wp) 
    709          zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpk-1,wp) 
    710          gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 
    711          gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 
     730         zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) 
     731         zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) 
     732         gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 
     733         gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 
    712734         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 
    713735      END DO 
     
    721743              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) ) 
    722744              ! 
    723               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) 
     745              e3w_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) 
    724746              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) 
    725747              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) 
     
    744766      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
    745767      !!---------------------------------------------------------------------- 
     768      USE utils, ONLY : wp 
    746769      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate 
    747770      REAL(wp)             ::   pf   ! sigma value 
    748771      !!---------------------------------------------------------------------- 
    749772      ! 
    750       pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpk-1) + rn_thetb )  )   & 
     773      pf =   (   TANH( rn_theta * ( -(pk-0.5) / REAL(jpk-1) + rn_thetb )  )   & 
    751774         &     - TANH( rn_thetb * rn_theta                                )  )   & 
    752775         & * (   COSH( rn_theta                           )                      & 
    753          &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              & 
    754          & / ( 2._wp * SINH( rn_theta ) ) 
     776         &     + COSH( rn_theta * ( 2. * rn_thetb - 1. ) )  )              & 
     777         & / ( 2. * SINH( rn_theta ) ) 
    755778      ! 
    756779   END FUNCTION fssig 
     
    768791      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
    769792      !!---------------------------------------------------------------------- 
     793      USE utils, ONLY : wp 
    770794      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate 
    771795      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient 
     
    774798      ! 
    775799      IF ( rn_theta == 0 ) then      ! uniform sigma 
    776          pf1 = - ( pk1 - 0.5_wp ) / REAL( jpk-1 ) 
     800         pf1 = - ( pk1 - 0.5 ) / REAL( jpk-1 ) 
    777801      ELSE                        ! stretched sigma 
    778          pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpk-1)) ) ) / SINH( rn_theta )              & 
    779             &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpk-1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  & 
    780             &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  ) 
     802         pf1 =   ( 1. - pbb ) * ( SINH( rn_theta*(-(pk1-0.5)/REAL(jpk-1)) ) ) / SINH( rn_theta )              & 
     803            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5)/REAL(jpk-1)) + 0.5) ) - TANH( 0.5 * rn_theta )  )  & 
     804            &        / ( 2. * TANH( 0.5 * rn_theta ) )  ) 
    781805      ENDIF 
    782806      ! 
     
    801825      !! Reference  :   Siddorn and Furner, in prep 
    802826      !!---------------------------------------------------------------------- 
     827      USE utils, ONLY : jpk,jk,wp 
    803828      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate 
    804829      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate 
     
    809834      REAL(wp)                ::   zn1,zn2        ! local variables 
    810835      REAL(wp)                ::   za,zb,zx       ! local variables 
    811       integer                 ::   jk 
    812836      !!---------------------------------------------------------------------- 
    813837      ! 
     
    816840      zn2  =  1. -  zn1 
    817841 
    818       za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp)  
    819       za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 
    820       za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 
     842      za1 = (rn_alpha+2.0)*zn1**(rn_alpha+1.0)-(rn_alpha+1.0)*zn1**(rn_alpha+2.0)  
     843      za2 = (rn_alpha+2.0)*zn2**(rn_alpha+1.0)-(rn_alpha+1.0)*zn2**(rn_alpha+2.0) 
     844      za3 = (zn2**3.0 - za2)/( zn1**3.0 - za1) 
    821845      
    822846      za = pzb - za3*(pzs-za1)-za2 
    823       za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 
    824       zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 
    825       zx = 1.0_wp-za/2.0_wp-zb 
     847      za = za/( zn2-0.5*(za2+zn2**2.0) - za3*(zn1-0.5*(za1+zn1**2.0) ) ) 
     848      zb = (pzs - za1 - za*( zn1-0.5*(za1+zn1**2.0 ) ) ) / (zn1**3.0 - za1) 
     849      zx = 1.0-za/2.0-zb 
    826850  
    827851      DO jk = 1, jpk 
    828         p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  & 
    829                     & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- & 
    830                     &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 
    831         p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 
     852        p_gamma(jk) = za*(pk1(jk)*(1.0-pk1(jk)/2.0))+zb*pk1(jk)**3.0 +  & 
     853                    & zx*( (rn_alpha+2.0)*pk1(jk)**(rn_alpha+1.0)- & 
     854                    &      (rn_alpha+1.0)*pk1(jk)**(rn_alpha+2.0) ) 
     855        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0-psmth) 
    832856      ENDDO  
    833857 
     
    835859   END FUNCTION fgamma 
    836860 
    837 END PROGRAM SCOORD_GEN 
  • branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90

    r5255 r5257  
    11MODULE utils 
    22 
    3    IMPLICIT NONE 
    43   USE netcdf 
    54 
     5   IMPLICIT NONE 
    66   PUBLIC             ! allows the acces to par_oce when dom_oce is used 
    77   !                  ! exception to coding rules... to be suppressed ??? 
    88 
    9    PUBLIC dom_oce_alloc 
    10    PUBLIC read_bathy 
     9!   PUBLIC dom_oce_alloc 
    1110 
     11   INTEGER, PARAMETER   :: dp=8 , sp=4, wp=dp 
    1212 
    1313   !! All coordinates 
     
    4545   INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers 
    4646   INTEGER  ::   ios                      ! Local integer output status for namelist read and allocation 
     47   INTEGER  ::   numnam                   ! File handle for namelist  
    4748   REAL(wp) ::   zrmax, ztaper   ! temporary scalars 
    4849   REAL(wp) ::   zrfact 
    4950   ! 
    50    REAL(wp), DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 
    51    REAL(wp), DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 
     51   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 
     52   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zenv, ztmp, zmsk, zri, zrj, zhbat 
     53 
     54   !Namelist variables 
     55   REAL(wp) :: rn_jpk, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax, rn_theta 
     56   REAL(wp) :: rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 
     57   LOGICAL :: ln_s_sh94, ln_s_sf12, ln_sigcrit 
    5258 
    5359   NAMELIST/namzgr_sco/rn_jpk, ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 
    5460                        rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 
    5561 
     62   CONTAINS 
    5663 
    5764   INTEGER FUNCTION dom_oce_alloc() 
    5865      !!---------------------------------------------------------------------- 
    59       INTEGER, DIMENSION(12) :: ierr 
     66      INTEGER, DIMENSION(4) :: ierr 
    6067      !!---------------------------------------------------------------------- 
    6168      ierr(:) = 0 
    6269      ! 
    6370      ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), & 
    64          &      zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj) ) 
     71         &      zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj), STAT=ierr(1) ) 
    6572         ! 
    66       ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) ,                         & 
    67          &      gdept_0 (jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) ,                         & 
    68          &      gdepw_0 (jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(4) ) 
     73      ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) ,                         & 
     74         &      gdept_0(jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) ,                         & 
     75         &      gdepw_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(2) ) 
    6976         ! 
    7077         ! 
     
    7481         &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     & 
    7582         &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     & 
    76          &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) ) 
     83         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(3) ) 
    7784 
    78       ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(9) ) 
    79       ! 
     85      ALLOCATE( mbathy(jpi,jpj) , STAT=ierr(4) ) 
     86     ! 
    8087      dom_oce_alloc = MAXVAL(ierr) 
    8188      ! 
    8289   END FUNCTION dom_oce_alloc 
    83  
     90    
    8491 
    8592   SUBROUTINE read_bathy() 
    8693     !! Read bathymetry from input netcdf file 
    87      INTEGER :: var_id 
     94     INTEGER :: var_id, ncin 
    8895 
    89      CALL check_nf90( nf90_open('bathy.nc', NF90_NOWRITE, ncin), 'Error opening mesh_mask file' ) 
     96     CALL check_nf90( nf90_open('bathy.nc', NF90_NOWRITE, ncin), 'Error opening bathy.nc file' ) 
    9097 
    9198     ! Find the size of the input bathymetry 
     
    96103      
    97104     ! Read the bathymetry variable from file 
    98      CALL check_nf90( nf90_inq_varid( ncin, 'bathymetry', tmask_id ), 'Cannot get variable ID for bathymetry') 
     105     CALL check_nf90( nf90_inq_varid( ncin, 'bathymetry', var_id ), 'Cannot get variable ID for bathymetry') 
    99106     CALL check_nf90( nf90_get_var( ncin, var_id, bathy, (/ 1,1 /), (/ jpi, jpj /) ) ) 
     107 
     108     CALL check_nf90( nf90_close(ncin), 'Error closing bathy.nc file' ) 
    100109 
    101110   END SUBROUTINE read_bathy 
     
    113122     CALL check_nf90( nf90_inquire_dimension(ncid,id_var,len=len)) 
    114123 
    115   END SUBROUTINE dimlen 
     124   END SUBROUTINE dimlen 
     125   
     126  
     127   SUBROUTINE write_coord_file() 
     128     ! Write out variables to the a netcdf coordinates file 
     129      
     130     INTEGER :: id_x, id_y, id_z 
     131     INTEGER :: ncout 
     132     INTEGER, DIMENSION(11) :: var_ids  !Array to contain all variable IDs 
    116133 
     134     !Create the file 
     135     CALL check_nf90( nf90_create('coord_zgr.nc', NF90_CLOBBER, ncout), 'Could not create output file') 
     136     ! 
     137     !Define dimensions 
     138     CALL check_nf90( nf90_def_dim(ncout, 'x', jpi, id_x) ) 
     139     CALL check_nf90( nf90_def_dim(ncout, 'y', jpj, id_y) ) 
     140     CALL check_nf90( nf90_def_dim(ncout, 'z', jpk, id_z) ) 
     141     ! 
     142     !Define variables 
     143     CALL check_nf90( nf90_def_var(ncout, 'gdept_0', nf90_double, (/id_x, id_y,id_x/), var_ids(1)) ) 
     144     CALL check_nf90( nf90_def_var(ncout, 'gdepw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(2)) ) 
     145     CALL check_nf90( nf90_def_var(ncout, 'gdep3w_0', nf90_double, (/id_x, id_y,id_x/), var_ids(3)) ) 
     146     CALL check_nf90( nf90_def_var(ncout, 'e3f_0', nf90_double, (/id_x, id_y,id_x/), var_ids(4)) ) 
     147     CALL check_nf90( nf90_def_var(ncout, 'e3t_0', nf90_double, (/id_x, id_y,id_x/), var_ids(5)) ) 
     148     CALL check_nf90( nf90_def_var(ncout, 'e3u_0', nf90_double, (/id_x, id_y,id_x/), var_ids(6)) ) 
     149     CALL check_nf90( nf90_def_var(ncout, 'e3v_0', nf90_double, (/id_x, id_y,id_x/), var_ids(7)) ) 
     150     CALL check_nf90( nf90_def_var(ncout, 'e3w_0', nf90_double, (/id_x, id_y,id_x/), var_ids(8)) ) 
     151     CALL check_nf90( nf90_def_var(ncout, 'e3uw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(9)) ) 
     152     CALL check_nf90( nf90_def_var(ncout, 'e3vw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(10)) ) 
     153     CALL check_nf90( nf90_def_var(ncout, 'mbathy', nf90_double, (/id_x, id_y,id_x/), var_ids(11)) ) 
     154      
     155     ! End define mode 
     156     CALL check_nf90( nf90_enddef(ncout) ) 
    117157 
    118     
    119    SUBROUTINE write_coord_file() 
     158     !Write variables to file 
     159     CALL check_nf90( nf90_put_var(ncout, var_ids(1), gdept_0) ) 
     160     CALL check_nf90( nf90_put_var(ncout, var_ids(2), gdepw_0) ) 
     161     CALL check_nf90( nf90_put_var(ncout, var_ids(3), gdep3w_0) ) 
     162     CALL check_nf90( nf90_put_var(ncout, var_ids(4), e3f_0) ) 
     163     CALL check_nf90( nf90_put_var(ncout, var_ids(5), e3t_0) ) 
     164     CALL check_nf90( nf90_put_var(ncout, var_ids(6), e3u_0) ) 
     165     CALL check_nf90( nf90_put_var(ncout, var_ids(7), e3v_0) ) 
     166     CALL check_nf90( nf90_put_var(ncout, var_ids(8), e3w_0) ) 
     167     CALL check_nf90( nf90_put_var(ncout, var_ids(9), e3uw_0) ) 
     168     CALL check_nf90( nf90_put_var(ncout, var_ids(10), e3vw_0) ) 
     169     CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy) ) 
     170      
     171     CALL check_nf90( nf90_close(ncout) ) 
    120172 
    121173   END SUBROUTINE write_coord_file 
     
    127179 
    128180      IF (istat /= nf90_noerr) THEN 
    129          WRITE(numerr,*) 'ERROR! : '//TRIM(nf90_strerror(istat)) 
    130          IF ( PRESENT(message) ) THEN ; WRITE(numerr,*) message ; ENDIF 
     181         WRITE(*,*) 'ERROR! : '//TRIM(nf90_strerror(istat)) 
     182         IF ( PRESENT(message) ) THEN ; WRITE(*,*) message ; ENDIF 
    131183         STOP 
    132184      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.