Changeset 1083


Ignore:
Timestamp:
2008-06-05T16:26:59+02:00 (13 years ago)
Author:
ctlod
Message:

trunk: small reorgansition of domzgr.F90 module, see ticket: #191

Location:
trunk/NEMO/OPA_SRC/DOM
Files:
1 deleted
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/DOM/domzgr.F90

    r900 r1083  
    1414   !!       zgr_zps      : z-coordinate with partial steps 
    1515   !!       zgr_sco      : s-coordinate 
     16   !!       fssig        : sigma coordinate non-dimensional function 
    1617   !!--------------------------------------------------------------------- 
    1718   !! * Modules used 
     
    729730   END SUBROUTINE zgr_zco 
    730731 
    731  
    732  
    733 #  include "domzgr_zps.h90" 
    734  
    735  
    736 #if ! defined key_zco 
     732#if defined key_zco 
    737733   !!---------------------------------------------------------------------- 
    738    !!   NOT 'key_zco' :                                    3D gdep. and e3. 
     734   !!   'key_zco' :                                              "pure" zco (gdep & e3 are 1D arrays) 
    739735   !!---------------------------------------------------------------------- 
     736   SUBROUTINE zgr_zps      ! Empty routine 
     737   END SUBROUTINE zgr_zps 
     738   SUBROUTINE zgr_sco      ! Empty routine 
     739   END SUBROUTINE zgr_sco 
     740 
     741#else 
     742   !!---------------------------------------------------------------------- 
     743   !!   Default option :                      zco, zps and/or sco available (gedp & e3 are 3D arrays) 
     744   !!---------------------------------------------------------------------- 
     745 
     746   SUBROUTINE zgr_zps 
     747      !!---------------------------------------------------------------------- 
     748      !!                  ***  ROUTINE zgr_zps  *** 
     749      !!                      
     750      !! ** Purpose :   the depth and vertical scale factor in partial step 
     751      !!      z-coordinate case 
     752      !! 
     753      !! ** Method  :   Partial steps : computes the 3D vertical scale factors 
     754      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with 
     755      !!      a partial step representation of bottom topography. 
     756      !! 
     757      !!        The reference depth of model levels is defined from an analytical 
     758      !!      function the derivative of which gives the reference vertical 
     759      !!      scale factors. 
     760      !!        From  depth and scale factors reference, we compute there new value 
     761      !!      with partial steps  on 3d arrays ( i, j, k ). 
     762      !! 
     763      !!              w-level: gdepw(i,j,k)  = fsdep(k) 
     764      !!                       e3w(i,j,k) = dk(fsdep)(k)     = fse3(i,j,k) 
     765      !!              t-level: gdept(i,j,k)  = fsdep(k+0.5) 
     766      !!                       e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5) 
     767      !! 
     768      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc), 
     769      !!      we find the mbathy index of the depth at each grid point. 
     770      !!      This leads us to three cases: 
     771      !! 
     772      !!              - bathy = 0 => mbathy = 0 
     773      !!              - 1 < mbathy < jpkm1     
     774      !!              - bathy > gdepw(jpk) => mbathy = jpkm1   
     775      !! 
     776      !!        Then, for each case, we find the new depth at t- and w- levels 
     777      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-  
     778      !!      and f-points. 
     779      !!  
     780      !!        This routine is given as an example, it must be modified 
     781      !!      following the user s desiderata. nevertheless, the output as 
     782      !!      well as the way to compute the model levels and scale factors 
     783      !!      must be respected in order to insure second order accuracy 
     784      !!      schemes. 
     785      !! 
     786      !!         c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives 
     787      !!         - - - - - - -   gdept, gdepw and e3. are positives 
     788      !!       
     789      !!  Reference : 
     790      !!     Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 
     791      !! 
     792      !!  History : 
     793      !!    8.5  !  02-09 (A. Bozec, G. Madec)  F90: Free form and module 
     794      !!    9.0  !  02-09 (A. de Miranda)  rigid-lid + islands 
     795      !!---------------------------------------------------------------------- 
     796      !! * Local declarations 
     797      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
     798      INTEGER  ::   & 
     799         ik, it            ! temporary integers 
     800       
     801      REAL(wp) ::   &   
     802         ze3tp, ze3wp,    &  ! Last ocean level thickness at T- and W-points 
     803         zdepwp,          &  ! Ajusted ocean depth to avoid too small e3t 
     804         zdepth,          &  !    "         " 
     805         zmax, zmin,      &  ! Maximum and minimum depth 
     806         zdiff               ! temporary scalar 
     807 
     808      REAL(wp), DIMENSION(jpi,jpj) ::   & 
     809         zprt  !    "           " 
     810 
     811      LOGICAL ::  ll_print                          ! Allow  control print for debugging 
     812 
     813      !!--------------------------------------------------------------------- 
     814      !! OPA8.5, LODYC-IPSL (2002) 
     815      !!--------------------------------------------------------------------- 
     816       
     817      ! Local variable for debugging 
     818      ll_print=.FALSE. 
     819!!!   ll_print=.TRUE. 
     820       
     821      ! Initialization of constant 
     822      zmax = gdepw_0(jpk) + e3t_0(jpk) 
     823      zmin = gdepw_0(4) 
     824       
     825      ! Ocean depth 
     826      IF(lwp .AND. ll_print) THEN  
     827         WRITE(numout,*) 
     828         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
     829         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
     830      ENDIF 
     831 
     832      IF(lwp) WRITE(numout,*) 
     833      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps' 
     834      IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
     835      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
     836 
     837 
     838      ! bathymetry in level (from bathy_meter) 
     839      ! =================== 
     840 
     841      ! initialize mbathy to the maximum ocean level available 
     842      mbathy(:,:) = jpkm1 
     843 
     844      ! storage of land and island's number (zera and negative values) in mbathy 
     845      DO jj = 1, jpj 
     846         DO ji= 1, jpi 
     847            IF( bathy(ji,jj) <= 0. )   mbathy(ji,jj) = INT( bathy(ji,jj) ) 
     848         END DO 
     849      END DO 
     850 
     851      ! bounded value of bathy 
     852      ! minimum depth == 3 levels 
     853      ! maximum depth == gdepw_0(jpk)+e3t_0(jpk)  
     854      ! i.e. the last ocean level thickness cannot exceed e3t_0(jpkm1)+e3t_0(jpk) 
     855      DO jj = 1, jpj 
     856         DO ji= 1, jpi 
     857            IF( bathy(ji,jj) <= 0. ) THEN 
     858               bathy(ji,jj) = 0.e0 
     859            ELSE 
     860               bathy(ji,jj) = MAX( bathy(ji,jj), zmin ) 
     861               bathy(ji,jj) = MIN( bathy(ji,jj), zmax ) 
     862            ENDIF 
     863         END DO 
     864      END DO 
     865 
     866      ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
     867      ! find the number of ocean levels such that the last level thickness 
     868      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where 
     869      ! e3t_0 is the reference level thickness 
     870      DO jk = jpkm1, 1, -1 
     871         zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat ) 
     872         DO jj = 1, jpj 
     873            DO ji = 1, jpi 
     874               IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth )   mbathy(ji,jj) = jk-1 
     875            END DO 
     876         END DO 
     877      END DO 
     878 
     879     ! Scale factors and depth at T- and W-points 
     880      
     881     ! intitialization to the reference z-coordinate 
     882     DO jk = 1, jpk 
     883        gdept(:,:,jk) = gdept_0(jk) 
     884        gdepw(:,:,jk) = gdepw_0(jk) 
     885        e3t  (:,:,jk) = e3t_0  (jk) 
     886        e3w  (:,:,jk) = e3w_0  (jk) 
     887     END DO 
     888     hdept(:,:) = gdept(:,:,2 ) 
     889     hdepw(:,:) = gdepw(:,:,3 )      
     890      
     891     !  
     892     DO jj = 1, jpj 
     893        DO ji = 1, jpi 
     894           ik = mbathy(ji,jj) 
     895           ! ocean point only 
     896           IF( ik > 0 ) THEN 
     897              ! max ocean level case 
     898              IF( ik == jpkm1 ) THEN 
     899                 zdepwp = bathy(ji,jj) 
     900                 ze3tp  = bathy(ji,jj) - gdepw_0(ik) 
     901                 ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) ) 
     902                 e3t(ji,jj,ik  ) = ze3tp 
     903                 e3t(ji,jj,ik+1) = ze3tp 
     904                 e3w(ji,jj,ik  ) = ze3wp 
     905                 e3w(ji,jj,ik+1) = ze3tp 
     906                 gdepw(ji,jj,ik+1) = zdepwp 
     907                 gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp 
     908                 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp 
     909                 ! standard case 
     910              ELSE 
     911!!alex 
     912                 IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN 
     913                    gdepw(ji,jj,ik+1) = bathy(ji,jj) 
     914                 ELSE 
     915!!alex ctl          write(*,*) 'zps',ji,jj,'bathy', bathy(ji,jj), 'depw ',gdepw_0(ik+1) 
     916                    gdepw(ji,jj,ik+1) = gdepw_0(ik+1) 
     917                 ENDIF 
     918!!Alex 
     919!!Alex           gdepw(ji,jj,ik+1) = bathy(ji,jj) 
     920!!bug gm  verifier les gdepw_0 
     921                 gdept(ji,jj,ik  ) =  gdepw_0(ik) + ( gdepw(ji,jj,ik+1) - gdepw_0(ik) )   & 
     922                    &                             * ((gdept_0(      ik  ) - gdepw_0(ik) )   & 
     923                    &                             / ( gdepw_0(      ik+1) - gdepw_0(ik) )) 
     924                 e3t(ji,jj,ik) = e3t_0(ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   &  
     925                    &                      / ( gdepw_0(      ik+1) - gdepw_0(ik) )  
     926                 e3w(ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) )   & 
     927                    &                * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) 
     928                 !       ... on ik+1 
     929                 e3w(ji,jj,ik+1) = e3t(ji,jj,ik) 
     930                 e3t(ji,jj,ik+1) = e3t(ji,jj,ik) 
     931                 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik) 
     932              ENDIF 
     933           ENDIF 
     934        END DO 
     935     END DO 
     936 
     937     it = 0 
     938     DO jj = 1, jpj 
     939        DO ji = 1, jpi 
     940           ik = mbathy(ji,jj) 
     941           ! ocean point only 
     942           IF( ik > 0 ) THEN 
     943              ! bathymetry output 
     944              hdept(ji,jj) = gdept(ji,jj,ik  ) 
     945              hdepw(ji,jj) = gdepw(ji,jj,ik+1) 
     946              e3tp (ji,jj) = e3t(ji,jj,ik  ) 
     947              e3wp (ji,jj) = e3w(ji,jj,ik  ) 
     948              ! test 
     949              zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  ) 
     950              IF( zdiff <= 0. .AND. lwp ) THEN  
     951                 it=it+1 
     952                 WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
     953                 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
     954                 WRITE(numout,*) ' gdept= ', gdept(ji,jj,ik), ' gdepw= ', gdepw(ji,jj,ik+1),   & 
     955                                 ' zdiff   = ', zdiff 
     956                 WRITE(numout,*) ' e3tp    = ', e3t(ji,jj,ik  ), ' e3wp    = ', e3w(ji,jj,ik  ) 
     957              ENDIF 
     958           ENDIF 
     959        END DO 
     960      END DO 
     961 
     962      ! Scale factors and depth at U-, V-, UW and VW-points 
     963 
     964      ! initialisation to z-scale factors 
     965      DO jk = 1, jpk 
     966         e3u (:,:,jk)  = e3t_0(jk) 
     967         e3v (:,:,jk)  = e3t_0(jk) 
     968         e3uw(:,:,jk)  = e3w_0(jk) 
     969         e3vw(:,:,jk)  = e3w_0(jk) 
     970      END DO 
     971 
     972     ! Computed as the minimum of neighbooring scale factors 
     973     DO jk = 1,jpk 
     974        DO jj = 1, jpjm1 
     975           DO ji = 1, fs_jpim1   ! vector opt. 
     976              e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk)) 
     977              e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk)) 
     978              e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) 
     979              e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) 
     980           END DO 
     981        END DO 
     982     END DO 
     983      
     984     ! Boundary conditions 
     985     CALL lbc_lnk( e3u , 'U', 1. )   ;   CALL lbc_lnk( e3uw, 'U', 1. ) 
     986     CALL lbc_lnk( e3v , 'V', 1. )   ;   CALL lbc_lnk( e3vw, 'V', 1. ) 
     987      
     988     ! set to z-scale factor if zero (i.e. along closed boundaries) 
     989     DO jk = 1, jpk 
     990        DO jj = 1, jpj 
     991           DO ji = 1, jpi 
     992              IF( e3u (ji,jj,jk) == 0.e0 ) e3u (ji,jj,jk) = e3t_0(jk) 
     993              IF( e3v (ji,jj,jk) == 0.e0 ) e3v (ji,jj,jk) = e3t_0(jk) 
     994              IF( e3uw(ji,jj,jk) == 0.e0 ) e3uw(ji,jj,jk) = e3w_0(jk) 
     995              IF( e3vw(ji,jj,jk) == 0.e0 ) e3vw(ji,jj,jk) = e3w_0(jk) 
     996           END DO 
     997        END DO 
     998     END DO 
     999      
     1000     ! Scale factor at F-point 
     1001      
     1002     ! initialisation to z-scale factors 
     1003     DO jk = 1, jpk 
     1004        e3f (:,:,jk) = e3t_0(jk) 
     1005     END DO 
     1006      
     1007     ! Computed as the minimum of neighbooring V-scale factors 
     1008     DO jk = 1, jpk 
     1009        DO jj = 1, jpjm1 
     1010           DO ji = 1, fs_jpim1   ! vector opt. 
     1011              e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) ) 
     1012           END DO 
     1013        END DO 
     1014     END DO 
     1015     ! Boundary conditions 
     1016     CALL lbc_lnk( e3f, 'F', 1. ) 
     1017      
     1018     ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1019     DO jk = 1, jpk 
     1020        DO jj = 1, jpj 
     1021           DO ji = 1, jpi 
     1022              IF( e3f(ji,jj,jk) == 0.e0 ) e3f(ji,jj,jk) = e3t_0(jk) 
     1023           END DO 
     1024        END DO 
     1025     END DO 
     1026!!bug ? gm:  must be a do loop with mj0,mj1 
     1027 
     1028     ! we duplicate factor scales for jj = 1 and jj = 2 
     1029     e3t(:,mj0(1),:) = e3t(:,mj0(2),:)  
     1030     e3w(:,mj0(1),:) = e3w(:,mj0(2),:)  
     1031     e3u(:,mj0(1),:) = e3u(:,mj0(2),:)  
     1032     e3v(:,mj0(1),:) = e3v(:,mj0(2),:)  
     1033     e3f(:,mj0(1),:) = e3f(:,mj0(2),:)  
     1034 
     1035 
     1036      
     1037     ! Compute gdep3w (vertical sum of e3w) 
     1038      
     1039     gdep3w   (:,:,1) = 0.5 * e3w (:,:,1) 
     1040     DO jk = 2, jpk 
     1041        gdep3w   (:,:,jk) = gdep3w   (:,:,jk-1) + e3w (:,:,jk)  
     1042     END DO 
     1043           
     1044     ! Control print 
     1045 9600 FORMAT(9x,' level   gdept    gdepw     e3t      e3w   ') 
     1046 9610 FORMAT(10x,i4,4f9.2) 
     1047      IF(lwp .AND. ll_print) THEN 
     1048         DO jj = 1,jpj 
     1049            DO ji = 1, jpi 
     1050               ik = MAX(mbathy(ji,jj),1) 
     1051               zprt(ji,jj) = e3t(ji,jj,ik) 
     1052            END DO 
     1053         END DO 
     1054         WRITE(numout,*) 
     1055         WRITE(numout,*) 'domzgr e3t(mbathy)' 
     1056         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1057         DO jj = 1,jpj 
     1058            DO ji = 1, jpi 
     1059               ik = MAX(mbathy(ji,jj),1) 
     1060               zprt(ji,jj) = e3w(ji,jj,ik) 
     1061            END DO 
     1062         END DO 
     1063         WRITE(numout,*) 
     1064         WRITE(numout,*) 'domzgr e3w(mbathy)' 
     1065         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1066         DO jj = 1,jpj 
     1067            DO ji = 1, jpi 
     1068               ik = MAX(mbathy(ji,jj),1) 
     1069               zprt(ji,jj) = e3u(ji,jj,ik) 
     1070            END DO 
     1071         END DO 
     1072 
     1073         WRITE(numout,*) 
     1074         WRITE(numout,*) 'domzgr e3u(mbathy)' 
     1075         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1076         DO jj = 1,jpj 
     1077            DO ji = 1, jpi 
     1078               ik = MAX(mbathy(ji,jj),1) 
     1079               zprt(ji,jj) = e3v(ji,jj,ik) 
     1080            END DO 
     1081         END DO 
     1082         WRITE(numout,*) 
     1083         WRITE(numout,*) 'domzgr e3v(mbathy)' 
     1084         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1085         DO jj = 1,jpj 
     1086            DO ji = 1, jpi 
     1087               ik = MAX(mbathy(ji,jj),1)  
     1088               zprt(ji,jj) = e3f(ji,jj,ik) 
     1089            END DO 
     1090         END DO 
     1091 
     1092         WRITE(numout,*) 
     1093         WRITE(numout,*) 'domzgr e3f(mbathy)' 
     1094         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1095         DO jj = 1,jpj 
     1096            DO ji = 1, jpi 
     1097               ik =  MAX(mbathy(ji,jj),1) 
     1098               zprt(ji,jj) = gdep3w(ji,jj,ik) 
     1099            END DO 
     1100         END DO 
     1101         WRITE(numout,*) 
     1102         WRITE(numout,*) 'domzgr gdep3w(mbathy)' 
     1103         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1104 
     1105      ENDIF   
     1106 
     1107 
     1108      DO jk = 1,jpk 
     1109         DO jj = 1, jpj 
     1110            DO ji = 1, jpi  
     1111               IF( e3w(ji,jj,jk) <= 0. .or. e3t(ji,jj,jk) <= 0. ) THEN 
     1112                  IF(lwp) THEN 
     1113                     WRITE(numout,*) ' e r r o r :         e3w or e3t =< 0 ' 
     1114                     WRITE(numout,*) ' =========           --------------- ' 
     1115                     WRITE(numout,*) 
     1116                  ENDIF 
     1117                  STOP 'domzgr.psteps' 
     1118               ENDIF 
     1119               IF( gdepw(ji,jj,jk) < 0. .or. gdept(ji,jj,jk) < 0. ) THEN 
     1120                  IF(lwp) THEN 
     1121                     WRITE(numout,*) ' e r r o r :      gdepw or gdept < 0 ' 
     1122                     WRITE(numout,*) ' =========        ------------------ ' 
     1123                     WRITE(numout,*) 
     1124                  ENDIF 
     1125                  STOP 'domzgr.psteps' 
     1126               ENDIF 
     1127            END DO 
     1128         END DO 
     1129      END DO   
     1130 
     1131   IF(lwp) THEN 
     1132      WRITE(numout,*) ' e3t lev 21 ' 
     1133      CALL prihre(e3t(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 
     1134      WRITE(numout,*) ' e3w lev 21  ' 
     1135      CALL prihre(e3w(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 
     1136      WRITE(numout,*) ' e3u lev 21  ' 
     1137      CALL prihre(e3u(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 
     1138      WRITE(numout,*) ' e3v lev 21  ' 
     1139      CALL prihre(e3v(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 
     1140      WRITE(numout,*) ' e3f lev 21  ' 
     1141      CALL prihre(e3f(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 
     1142      WRITE(numout,*) ' e3t lev 22 ' 
     1143      CALL prihre(e3t(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 
     1144      WRITE(numout,*) ' e3w lev 22  ' 
     1145      CALL prihre(e3w(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 
     1146      WRITE(numout,*) ' e3u lev 22  ' 
     1147      CALL prihre(e3u(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 
     1148      WRITE(numout,*) ' e3v lev 22  ' 
     1149      CALL prihre(e3v(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 
     1150      WRITE(numout,*) ' e3f lev 22  ' 
     1151      CALL prihre(e3f(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 
     1152   ENDIF 
     1153 
     1154      ! =========== 
     1155      ! Zoom domain  
     1156      ! =========== 
     1157 
     1158      IF( lzoom )   CALL zgr_bat_zoom 
     1159 
     1160      ! ================ 
     1161      ! Bathymetry check 
     1162      ! ================ 
     1163 
     1164      CALL zgr_bat_ctl 
     1165 
     1166   END SUBROUTINE zgr_zps 
     1167 
     1168 
     1169   FUNCTION fssig( pk ) RESULT( pf ) 
     1170      !!---------------------------------------------------------------------- 
     1171      !!                 ***  ROUTINE eos_init  *** 
     1172      !!        
     1173      !! ** Purpose :   provide the analytical function in s-coordinate 
     1174      !!           
     1175      !! ** Method  :   the function provide the non-dimensional position of 
     1176      !!                T and W (i.e. between 0 and 1) 
     1177      !!                T-points at integer values (between 1 and jpk) 
     1178      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
     1179      !! 
     1180      !! Reference  :   ??? 
     1181      !!---------------------------------------------------------------------- 
     1182      REAL(wp), INTENT(in   ) ::   pk   ! continuous "k" coordinate 
     1183      REAL(wp)                ::   pf   ! sigma value 
     1184      !!---------------------------------------------------------------------- 
     1185      ! 
     1186      pf =   (   TANH( theta * ( -(pk-0.5) / REAL(jpkm1) + thetb )  )      & 
     1187         &     - TANH( thetb * theta                                )  )   & 
     1188         & * (   COSH( theta                           )                   & 
     1189         &     + COSH( theta * ( 2.e0 * thetb - 1.e0 ) )  )                & 
     1190         & / ( 2.e0 * SINH( theta ) ) 
     1191      ! 
     1192   END FUNCTION fssig 
     1193 
    7401194 
    7411195   SUBROUTINE zgr_sco 
     
    7581212      !!            hbatf = mi( mj( hbatt ) ) 
    7591213      !!          - Compute gsigt, gsigw, esigt, esigw from an analytical 
    760       !!         function and its derivative given as statement function. 
    761       !!            gsigt(k) = fssig (k+0.5) 
    762       !!            gsigw(k) = fssig (k    ) 
    763       !!            esigt(k) = fsdsig(k+0.5) 
    764       !!            esigw(k) = fsdsig(k    ) 
     1214      !!         function and its derivative given as function. 
     1215      !!            gsigt(k) = fssig (k    ) 
     1216      !!            gsigw(k) = fssig (k-0.5) 
     1217      !!            esigt(k) = fsdsig(k    ) 
     1218      !!            esigw(k) = fsdsig(k-0.5) 
    7651219      !!      This routine is given as an example, it must be modified 
    7661220      !!      following the user s desiderata. nevertheless, the output as 
     
    7811235      !! * Local declarations 
    7821236      INTEGER  ::   ji, jj, jk, jl 
    783       REAL(wp) ::   fssig, fsdsig, pfkk 
    784  
    7851237      INTEGER  ::   iip1, ijp1, iim1, ijm1 
    786       REAL(wp) ::   & 
    787          fssigt, fssigw, zcoeft, zcoefw,   & 
    788          zrmax, ztaper 
     1238      REAL(wp) ::   zrmax, ztaper 
    7891239 
    7901240      REAL(wp), DIMENSION(jpi,jpj) ::   & 
     
    7921242 
    7931243      NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max 
    794       !!---------------------------------------------------------------------- 
    795  
    796       ! a. Sigma stretching coefficient 
    797        fssigt(pfkk) = (tanh(theta*(-(pfkk-0.5)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) & 
    798                     * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta)) 
    799        fssigw(pfkk) = (tanh(theta*(-(pfkk-1.0)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) & 
    800                     * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta)) 
    801  
    802       ! b. Vertical derivative of sigma stretching coefficient 
    803   
    804 !!bug fsdsig(pfkk)= put here the analytical derivative of fssigt and w ... 
    8051244      !!---------------------------------------------------------------------- 
    8061245 
     
    10431482      ! 2.1 Coefficients for model level depth at w- and t-levels 
    10441483 
    1045 !!bug gm : change the def of statment function.... 
    10461484      DO jk = 1, jpk 
    1047         gsigw(jk) = -fssigw(FLOAT(jk)) 
    1048         gsigt(jk) = -fssigt(FLOAT(jk)) 
     1485        gsigw(jk) = -fssig( FLOAT(jk)-0.5 ) 
     1486        gsigt(jk) = -fssig( FLOAT(jk)     ) 
    10491487      END DO 
    10501488 
    10511489      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk) 
    1052  
    1053 !!org DO jk = 1, jpk 
    1054 !!org    gsigw(jk) = -fssig ( FLOAT(jk)    ) 
    1055 !!org    gsigt(jk) = -fssig ( FLOAT(jk)+0.5) 
    1056 !!org END DO 
    10571490 
    10581491      ! 2.2 Coefficients for vertical scale factors at w-, t- levels 
     
    10681501        esigw(1)=esigw(2) 
    10691502        esigt(jpk)=esigt(jpkm1) 
    1070  
    1071 !!org DO jk = 1, jpk 
    1072 !!org    esigw(jk)=fsdsig( FLOAT(jk)) 
    1073 !!org    esigt(jk)=fsdsig( FLOAT(jk)+0.5) 
    1074 !!org END DO 
    10751503 
    10761504      ! 2.3 Coefficients for vertical depth as the sum of e3w scale factors 
     
    12051633   END SUBROUTINE zgr_sco 
    12061634 
    1207 #else 
    1208    !!---------------------------------------------------------------------- 
    1209    !!   Default option :                                      Empty routine 
    1210    !!---------------------------------------------------------------------- 
    1211    SUBROUTINE zgr_sco 
    1212    END SUBROUTINE zgr_sco 
    12131635#endif 
    12141636 
Note: See TracChangeset for help on using the changeset viewer.