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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11258 for NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src – NEMO

Ignore:
Timestamp:
2019-07-11T18:21:02+02:00 (5 years ago)
Author:
smasson
Message:

dev_r10984_HPC-13 : minor bugfix and cleaning, see #2285

Location:
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn2d.F90

    r11224 r11258  
    199199      ELSE                ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd) 
    200200      END IF 
    201       sshdta(:,:) = 0.0 
     201      ! 
    202202      DO jb = ibeg, iend 
    203203         ii = idx%nbi(jb,igrd) 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90

    r11223 r11258  
    116116         ! Open boundaries definition (arrays and masks) 
    117117         CALL bdy_def 
     118         IF( ln_meshmask )   CALL bdy_meshwri() 
    118119         ! 
    119120         ! Open boundaries initialisation of external data arrays 
     
    146147      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, ir, iseg     ! dummy loop indices 
    147148      INTEGER  ::   icount, icountr, icountr0, ibr_max     ! local integers 
    148       INTEGER  ::   ilen1, ibm1                            !   -       - 
     149      INTEGER  ::   ilen1                                  !   -       - 
    149150      INTEGER  ::   iwe, ies, iso, ino, inum, id_dummy     !   -       - 
    150151      INTEGER  ::   jpbdta                                 !   -       - 
     
    153154      INTEGER  ::   iibe, ijbe, iibi, ijbi                 !   -       - 
    154155      INTEGER  ::   flagu, flagv                           ! short cuts 
    155       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zz_read      ! work space for 2D global boundary data 
    156       REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask   ! pointer to 2D mask fields 
    157       INTEGER, DIMENSION (4)                  ::   kdimsz 
    158       INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays  
    159       INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta 
    160       INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points 
    161       CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid 
    162       REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                  ! temporary fmask array excluding coastal boundary condition (shlat) 
    163       REAL(wp), DIMENSION(jpi,jpj) ::   ztmask, zumask, zvmask  ! temporary u/v mask array 
    164       !! 
    165       INTEGER                              ::   nbdyind, nbdybeg, nbdyend 
     156      INTEGER  ::   nbdyind, nbdybeg, nbdyend 
     157      INTEGER              , DIMENSION(4)             ::   kdimsz 
     158      INTEGER              , DIMENSION(jpbgrd,jp_bdy) ::   nblendta          ! Length of index arrays  
     159      INTEGER,  ALLOCATABLE, DIMENSION(:,:,:)         ::   nbidta, nbjdta    ! Index arrays: i and j indices of bdy dta 
     160      INTEGER,  ALLOCATABLE, DIMENSION(:,:,:)         ::   nbrdta            ! Discrete distance from rim points 
     161      CHARACTER(LEN=1)     , DIMENSION(jpbgrd)        ::   cgrid 
     162      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zz_read                 ! work space for 2D global boundary data 
     163      REAL(wp), POINTER    , DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields 
     164      REAL(wp)             , DIMENSION(jpi,jpj) ::   zfmask   ! temporary fmask array excluding coastal boundary condition (shlat) 
     165      REAL(wp)             , DIMENSION(jpi,jpj) ::   ztmask, zumask, zvmask  ! temporary u/v mask array 
    166166      !!---------------------------------------------------------------------- 
    167167      ! 
     
    483483         END DO 
    484484      END DO 
    485  
     485      ! 
    486486      ! Find lenght of boundaries and rim on local mpi domain 
    487487      !------------------------------------------------------ 
     
    502502            DO ib = 1, nblendta(igrd,ib_bdy) 
    503503               ! check that data is in correct order in file 
    504                ibm1 = MAX(1,ib-1) 
    505                IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc... 
    506                   IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 
     504               IF( ib > 1 ) THEN 
     505                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ib-1,igrd,ib_bdy) ) THEN 
    507506                     CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', & 
    508507                        &        ' in order of distance from edge nbr A utility for re-ordering ', & 
     
    636635      ! For the flagu/flagv calculation below we require a version of fmask without 
    637636      ! the land boundary condition (shlat) included: 
    638       DO ij = 2, jpjm1 
    639          DO ii = 2, jpim1 
    640             zfmask(ii,ij) = ztmask(ii,ij  ) * ztmask(ii+1,ij  )   & 
    641                &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 
     637      DO ij = 1, jpjm1 
     638         DO ii = 1, jpim1 
     639            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   & 
     640               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 
    642641         END DO 
    643642      END DO 
     
    678677      ! handle rim1, do as if rim 0 was land 
    679678      ! ------------------------------------ 
    680  
     679       
     680      ! z[tuv]mask are now set to zero on rim 0 points: 
    681681      DO ib_bdy = 1, nb_bdy 
    682682         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0 
     
    692692 
    693693      ! Recompute zfmask 
    694       DO ij = 2, jpjm1 
    695          DO ii = 2, jpim1 
    696             zfmask(ii,ij) = ztmask(ii,ij  ) * ztmask(ii+1,ij  )   & 
    697                &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 
     694      DO ij = 1, jpjm1 
     695         DO ii = 1, jpim1 
     696            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   & 
     697               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 
    698698         END DO 
    699699      END DO 
     
    714714 
    715715      CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. )   ! compute flagu, flagv, ntreat on rim 1 
    716  
    717716      ! 
    718717      ! Check which boundaries might need communication 
     
    727726         DO ib_bdy = 1, nb_bdy 
    728727            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     728               IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE 
    729729               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    730730               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    731                IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN   ;   ir = 0 
    732                ELSE                                                 ;   ir = 1 
    733                END IF 
     731               ir = idx_bdy(ib_bdy)%nbr(ib,igrd) 
    734732               flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd)) 
    735733               flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd)) 
     
    784782 
    785783      DO ib_bdy = 1,nb_bdy 
    786          IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. & 
     784         IF(  cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. & 
    787785            & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. & 
    788786            & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN 
     
    793791                  IF(  mig(ii) > 2 .AND. mig(ii) < jpiglo-2 .AND. mjg(ij) > 2 .AND. mjg(ij) < jpjglo-2  ) THEN 
    794792                     WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain' 
    795                      WRITE(ctmp2,*) ' ========== ' 
    796                      CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     793                     CALL ctl_stop( ctmp1 ) 
    797794                  END IF 
    798795               END DO 
     
    800797         END IF 
    801798      END DO 
    802  
    803       ! 
    804       ! Tidy up 
    805       !-------- 
     799      ! 
    806800      DEALLOCATE( nbidta, nbjdta, nbrdta ) 
    807801      ! 
     
    809803 
    810804 
    811    SUBROUTINE bdy_rim_treat( zumask, zvmask, zfmask, lrim0 ) 
     805   SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 ) 
    812806      !!---------------------------------------------------------------------- 
    813807      !!                 ***  ROUTINE bdy_rim_treat  *** 
     
    826820      !!                - and look at the ocean neighbours to compute ntreat 
    827821      !!---------------------------------------------------------------------- 
    828       REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: zfmask   ! temporary fmask excluding coastal boundary condition (shlat) 
    829       REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: zumask, zvmask   ! temporary t/u/v mask array 
     822      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pfmask   ! temporary fmask excluding coastal boundary condition (shlat) 
     823      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pumask, pvmask   ! temporary t/u/v mask array 
    830824      LOGICAL                             , INTENT (in   ) :: lrim0    ! .true. -> rim 0   .false. -> rim 1 
    831825      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, icount       ! dummy loop indices 
     
    833827      INTEGER  ::   ibeg, iend                             ! local integer 
    834828      LOGICAL  ::   llnon, llson, llean, llwen             ! local logicals indicating the presence of a ocean neighbour 
    835       REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask   ! pointer to 2D mask fields 
     829      REAL(wp), POINTER, DIMENSION(:,:)       ::   zmask   ! pointer to 2D mask fields 
    836830      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
    837831      CHARACTER(LEN=1), DIMENSION(jpbgrd)     ::   cgrid 
     
    849843         DO igrd = 1, jpbgrd  
    850844            SELECT CASE( igrd ) 
    851                CASE( 1 )   ;   pmask => zumask     ;   i_offset = 0 
    852                CASE( 2 )   ;   pmask => bdytmask   ;   i_offset = 1 
    853                CASE( 3 )   ;   pmask => zfmask     ;   i_offset = 0 
     845               CASE( 1 )   ;   zmask => pumask     ;   i_offset = 0 
     846               CASE( 2 )   ;   zmask => bdytmask   ;   i_offset = 1 
     847               CASE( 3 )   ;   zmask => pfmask     ;   i_offset = 0 
    854848            END SELECT  
    855849            icount = 0 
    856             ztmp(:,:) = 0._wp 
     850            ztmp(:,:) = -999._wp 
    857851            IF( lrim0 ) THEN   ! extent of rim 0 
    858852               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 
     
    864858               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    865859               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE 
    866                zwfl = pmask(ii+i_offset-1,ij) 
    867                zefl = pmask(ii+i_offset  ,ij) 
     860               zwfl = zmask(ii+i_offset-1,ij) 
     861               zefl = zmask(ii+i_offset  ,ij) 
    868862               ! This error check only works if you are using the bdyXmask arrays 
    869863               IF( i_offset == 1 .and. zefl + zwfl == 2. ) THEN 
     
    875869            END DO 
    876870            IF( icount /= 0 ) THEN 
    877                WRITE(ctmp1,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   & 
     871               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   & 
    878872                  ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 
    879                WRITE(ctmp2,*) ' ========== ' 
    880                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     873               CALL ctl_stop( ctmp1 ) 
    881874            ENDIF  
    882             CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )  
     875            SELECT CASE( igrd ) 
     876               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )  
     877               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. )  
     878               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. )  
     879            END SELECT  
    883880            DO ib = ibeg, iend 
    884881               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     
    894891         DO igrd = 1, jpbgrd  
    895892            SELECT CASE( igrd ) 
    896                CASE( 1 )   ;   pmask => zvmask     ;   j_offset = 0 
    897                CASE( 2 )   ;   pmask => zfmask     ;   j_offset = 0 
    898                CASE( 3 )   ;   pmask => bdytmask   ;   j_offset = 1 
     893               CASE( 1 )   ;   zmask => pvmask     ;   j_offset = 0 
     894               CASE( 2 )   ;   zmask => pfmask     ;   j_offset = 0 
     895               CASE( 3 )   ;   zmask => bdytmask   ;   j_offset = 1 
    899896            END SELECT  
    900897            icount = 0 
    901             ztmp(:,:) = 0._wp 
     898            ztmp(:,:) = -999._wp 
    902899            IF( lrim0 ) THEN   ! extent of rim 0 
    903900               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 
     
    909906               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    910907               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE 
    911                zsfl = pmask(ii,ij+j_offset-1) 
    912                znfl = pmask(ii,ij+j_offset  ) 
     908               zsfl = zmask(ii,ij+j_offset-1) 
     909               znfl = zmask(ii,ij+j_offset  ) 
    913910               ! This error check only works if you are using the bdyXmask arrays 
    914911               IF( j_offset == 1 .and. znfl + zsfl == 2. ) THEN 
     
    920917            END DO 
    921918            IF( icount /= 0 ) THEN 
    922                WRITE(ctmp1,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   & 
     919               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   & 
    923920                  ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 
    924                WRITE(ctmp2,*) ' ========== ' 
    925                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     921               CALL ctl_stop( ctmp1 ) 
    926922            ENDIF 
    927             CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )  
     923            SELECT CASE( igrd ) 
     924               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )  
     925               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. )  
     926               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. )  
     927            END SELECT  
    928928            DO ib = ibeg, iend 
    929929               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     
    938938         DO igrd = 1, jpbgrd 
    939939            SELECT CASE( igrd ) 
    940                CASE( 1 )   ;   pmask => bdytmask  
    941                CASE( 2 )   ;   pmask => bdyumask  
    942                CASE( 3 )   ;   pmask => bdyvmask  
     940               CASE( 1 )   ;   zmask => bdytmask  
     941               CASE( 2 )   ;   zmask => bdyumask  
     942               CASE( 3 )   ;   zmask => bdyvmask  
    943943            END SELECT 
    944             ztmp(:,:) = 0._wp 
     944            ztmp(:,:) = -999._wp 
    945945            IF( lrim0 ) THEN   ! extent of rim 0 
    946946               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 
     
    952952               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    953953               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE 
    954                llnon = pmask(ii  ,ij+1) == 1.   
    955                llson = pmask(ii  ,ij-1) == 1.  
    956                llean = pmask(ii+1,ij  ) == 1.  
    957                llwen = pmask(ii-1,ij  ) == 1.  
     954               llnon = zmask(ii  ,ij+1) == 1.   
     955               llson = zmask(ii  ,ij-1) == 1.  
     956               llean = zmask(ii+1,ij  ) == 1.  
     957               llwen = zmask(ii-1,ij  ) == 1.  
    958958               inn  = COUNT( (/ llnon, llson, llean, llwen /) ) 
    959959               IF( inn == 0 ) THEN   ! no neighbours -> interior of a corner  or  cluster of rim points 
     
    961961                  !  1 |   o      !  2  o   |    !  3 | x        !  4     x |   !      |   |   -> error 
    962962                  !    |_x_ _     !    _ _x_|    !    |   o      !      o   |   !      |x_x| 
    963                   IF(     pmask(ii+1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 1. 
    964                   ELSEIF( pmask(ii-1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 2. 
    965                   ELSEIF( pmask(ii+1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 3. 
    966                   ELSEIF( pmask(ii-1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 4. 
    967                   ELSE 
    968                      WRITE(ctmp1,*) ' E R R O R : Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   & 
     963                  IF(     zmask(ii+1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 1. 
     964                  ELSEIF( zmask(ii-1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 2. 
     965                  ELSEIF( zmask(ii+1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 3. 
     966                  ELSEIF( zmask(ii-1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 4. 
     967                  ELSE                                    ;   ztmp(ii,ij) = -1. 
     968                     WRITE(ctmp1,*) 'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   & 
    969969                       ' on boundary set ', ib_bdy, ' has no free ocean neighbour' 
    970970                     IF( lrim0 ) THEN 
     
    973973                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.' 
    974974                     END IF 
    975                      WRITE(ctmp3,*) ' ========== ' 
    976                      CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) 
     975                     CALL ctl_warn( ctmp1, ctmp2 ) 
    977976                  END IF 
    978977               END IF 
     
    10051004               END IF 
    10061005               IF( inn == 4 ) THEN 
    1007                   WRITE(ctmp1,*)  ' E R R O R : Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   & 
     1006                  WRITE(ctmp1,*)  'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   & 
    10081007                       ' on boundary set ', ib_bdy, ' have 4 neighbours' 
    1009                   WRITE(ctmp2,*) ' ========== ' 
    1010                   CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1008                  CALL ctl_stop( ctmp1 ) 
    10111009               END IF 
    10121010            END DO 
    1013             CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
     1011            SELECT CASE( igrd ) 
     1012               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )  
     1013               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. )  
     1014               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. )  
     1015            END SELECT  
    10141016            DO ib = ibeg, iend 
    10151017               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     
    10311033      !! 
    10321034      !! ** Method  :  use itreat input to select a case 
    1033       !!               N.B. ntreat is defined for all bdy points in routine bdy_segs 
     1035      !!               N.B. ntreat is defined for all bdy points in routine bdy_rim_treat 
    10341036      !! 
    10351037      !!---------------------------------------------------------------------- 
     
    12431245                     icorns(ib2,1) = npckgw(ib1) 
    12441246                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN 
    1245                      WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1247                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 
    12461248                        &                                     jpisft(ib2), jpjwft(ib1) 
    1247                      WRITE(ctmp2,*) ' ==========  Not allowed yet' 
    1248                      WRITE(ctmp3,*) '             Crossing problem with West segment: ',npckgw(ib1), &  
    1249                         &                                        ' and South segment: ',npckgs(ib2) 
    1250                      CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) 
     1249                     WRITE(ctmp2,*) ' Not allowed yet' 
     1250                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &  
     1251                        &                            ' and South segment: ',npckgs(ib2) 
     1252                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 
    12511253                  ELSE 
    1252                      WRITE(ctmp1,*) ' E R R O R : Check South and West Open boundary indices' 
    1253                      WRITE(ctmp2,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1) , & 
    1254                         &                                         ' and South segment: ',npckgs(ib2) 
    1255                      CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1254                     WRITE(ctmp1,*) ' Check South and West Open boundary indices' 
     1255                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , & 
     1256                        &                            ' and South segment: ',npckgs(ib2) 
     1257                     CALL ctl_stop( ctmp1, ctmp2 ) 
    12561258                  END IF 
    12571259               END IF 
     
    12751277                     icorns(ib2,2) = npckge(ib1) 
    12761278                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN 
    1277                      WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1279                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 
    12781280                        &                                     jpisdt(ib2), jpjeft(ib1) 
    1279                      WRITE(ctmp2,*) ' ==========  Not allowed yet' 
    1280                      WRITE(ctmp3,*) '             Crossing problem with East segment: ',npckge(ib1), & 
    1281                         &                                        ' and South segment: ',npckgs(ib2) 
    1282                      CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) 
     1281                     WRITE(ctmp2,*) ' Not allowed yet' 
     1282                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), & 
     1283                        &                            ' and South segment: ',npckgs(ib2) 
     1284                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 
    12831285                  ELSE 
    1284                      WRITE(ctmp1,*) ' E R R O R : Check South and East Open boundary indices' 
    1285                      WRITE(ctmp2,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), & 
    1286                      &                                           ' and South segment: ',npckgs(ib2) 
    1287                      CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1286                     WRITE(ctmp1,*) ' Check South and East Open boundary indices' 
     1287                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), & 
     1288                     &                               ' and South segment: ',npckgs(ib2) 
     1289                     CALL ctl_stop( ctmp1, ctmp2 ) 
    12881290                  END IF 
    12891291               END IF 
     
    13071309                     icornn(ib2,1) = npckgw(ib1) 
    13081310                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN 
    1309                      WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1311                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 
    13101312                     &                                     jpinft(ib2), jpjwdt(ib1) 
    1311                      WRITE(ctmp2,*) ' ==========  Not allowed yet' 
    1312                      WRITE(ctmp3,*) '             Crossing problem with West segment: ',npckgw(ib1), & 
    1313                      &                                                    ' and North segment: ',npckgn(ib2) 
    1314                      CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) 
     1313                     WRITE(ctmp2,*) ' Not allowed yet' 
     1314                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 
     1315                     &                               ' and North segment: ',npckgn(ib2) 
     1316                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 
    13151317                  ELSE 
    1316                      WRITE(ctmp1,*) ' E R R O R : Check North and West Open boundary indices' 
    1317                      WRITE(ctmp2,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1), & 
    1318                      &                                                    ' and North segment: ',npckgn(ib2) 
    1319                      CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1318                     WRITE(ctmp1,*) ' Check North and West Open boundary indices' 
     1319                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), & 
     1320                     &                               ' and North segment: ',npckgn(ib2) 
     1321                     CALL ctl_stop( ctmp1, ctmp2 ) 
    13201322                  END IF 
    13211323               END IF 
     
    13391341                     icornn(ib2,2) = npckge(ib1) 
    13401342                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN 
    1341                      WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1343                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 
    13421344                     &                                     jpindt(ib2), jpjedt(ib1) 
    1343                      WRITE(ctmp2,*) ' ==========  Not allowed yet' 
    1344                      WRITE(ctmp3,*) '             Crossing problem with East segment: ',npckge(ib1), & 
    1345                      &                                           ' and North segment: ',npckgn(ib2) 
    1346                      CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) 
     1345                     WRITE(ctmp2,*) ' Not allowed yet' 
     1346                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), & 
     1347                     &                               ' and North segment: ',npckgn(ib2) 
     1348                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 
    13471349                  ELSE 
    1348                      WRITE(ctmp1,*) ' E R R O R : Check North and East Open boundary indices' 
    1349                      WRITE(ctmp2,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), & 
    1350                      &                                           ' and North segment: ',npckgn(ib2) 
    1351                      CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1350                     WRITE(ctmp1,*) ' Check North and East Open boundary indices' 
     1351                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), & 
     1352                     &                               ' and North segment: ',npckgn(ib2) 
     1353                     CALL ctl_stop( ctmp1, ctmp2 ) 
    13521354                  END IF 
    13531355               END IF 
     
    13751377         IF (ztestmask(1)==1) THEN  
    13761378            IF (icornw(ib,1)==0) THEN 
    1377                WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib) 
    1378                WRITE(ctmp2,*) ' ==========  does not start on land or on a corner'                                                   
    1379                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1379               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib) 
     1380               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 
    13801381            ELSE 
    13811382               ! This is a corner 
     
    13871388         IF (ztestmask(2)==1) THEN 
    13881389            IF (icornw(ib,2)==0) THEN 
    1389                WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib) 
    1390                WRITE(ctmp2,*) ' ==========  does not end on land or on a corner'                                                   
    1391                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1390               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib) 
     1391               CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' ) 
    13921392            ELSE 
    13931393               ! This is a corner 
     
    14151415         IF (ztestmask(1)==1) THEN 
    14161416            IF (icorne(ib,1)==0) THEN 
    1417                WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib) 
    1418                WRITE(ctmp2,*) ' ==========  does not start on land or on a corner'                                                   
    1419                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1417               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib) 
     1418               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 
    14201419            ELSE 
    14211420               ! This is a corner 
     
    14271426         IF (ztestmask(2)==1) THEN 
    14281427            IF (icorne(ib,2)==0) THEN 
    1429                WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib) 
    1430                WRITE(ctmp2,*) ' ==========  does not end on land or on a corner'                                                   
    1431                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1428               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib) 
     1429               CALL ctl_stop( ctmp1, ' does not end on land or on a corner' ) 
    14321430            ELSE 
    14331431               ! This is a corner 
     
    14541452 
    14551453         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN 
    1456             WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib) 
    1457             WRITE(ctmp2,*) ' ==========  does not start on land or on a corner'                                                   
    1458             CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1454            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib) 
     1455            CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 
    14591456         ENDIF 
    14601457         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN 
    1461             WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib) 
    1462             WRITE(ctmp2,*) ' ==========  does not end on land or on a corner'                                                   
    1463             CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1458            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib) 
     1459            CALL ctl_stop( ctmp1, ' does not end on land or on a corner' ) 
    14641460         ENDIF 
    14651461      END DO 
     
    14801476 
    14811477         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN 
    1482             WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib) 
    1483             WRITE(ctmp2,*) ' ==========  does not start on land'                                                   
    1484             CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1478            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib) 
     1479            CALL ctl_stop( ctmp1, ' does not start on land' ) 
    14851480         ENDIF 
    14861481         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN 
    1487             WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib) 
    1488             WRITE(ctmp2,*) ' ==========  does not end on land'                                                   
    1489             CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1482            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib) 
     1483            CALL ctl_stop( ctmp1, ' does not end on land' ) 
    14901484         ENDIF 
    14911485      END DO 
     
    17271721      ! 
    17281722      IF( itest>0 ) THEN 
    1729          WRITE(ctmp1,*) ' E R R O R : Segments ', ib1, 'and ', ib2 
    1730          WRITE(ctmp2,*) ' ==========  have different open bdy schemes'                                                   
    1731          CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1723         WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2 
     1724         CALL ctl_stop( ctmp1, ' have different open bdy schemes' ) 
    17321725      ENDIF 
    17331726      ! 
    17341727   END SUBROUTINE bdy_ctl_corn 
    17351728 
     1729 
     1730   SUBROUTINE bdy_meshwri() 
     1731      !!---------------------------------------------------------------------- 
     1732      !!                 ***  ROUTINE bdy_meshwri  *** 
     1733      !!          
     1734      !! ** Purpose :   write netcdf file with nbr, flagu, flagv, ntreat for T, U  
     1735      !!                and V points in 2D arrays for easier visualisation/control 
     1736      !! 
     1737      !! ** Method  :   use iom_rstput as in domwri.F 
     1738      !!----------------------------------------------------------------------       
     1739      INTEGER  ::   ib_bdy, ii, ij, igrd, ib     ! dummy loop indices 
     1740      INTEGER  ::   inum                                   !   -       - 
     1741      REAL(wp), POINTER, DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields 
     1742      REAL(wp)         , DIMENSION(jpi,jpj) ::   ztmp 
     1743      CHARACTER(LEN=1) , DIMENSION(jpbgrd)  ::   cgrid 
     1744      !!----------------------------------------------------------------------       
     1745      cgrid = (/'t','u','v'/) 
     1746      CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. ) 
     1747      DO igrd = 1, jpbgrd 
     1748         SELECT CASE( igrd ) 
     1749         CASE( 1 )   ;   zmask => tmask(:,:,1) 
     1750         CASE( 2 )   ;   zmask => umask(:,:,1) 
     1751         CASE( 3 )   ;   zmask => vmask(:,:,1) 
     1752         END SELECT 
     1753         ztmp(:,:) = zmask(:,:) 
     1754         DO ib_bdy = 1, nb_bdy 
     1755            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)      ! nbr deined for all rims 
     1756               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     1757               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     1758               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10. 
     1759               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 
     1760            END DO 
     1761         END DO 
     1762         CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 
     1763         ztmp(:,:) = zmask(:,:) 
     1764         DO ib_bdy = 1, nb_bdy 
     1765            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagu defined only for rims 0 and 1 
     1766               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     1767               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     1768               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10. 
     1769               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 
     1770            END DO 
     1771         END DO 
     1772         CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 
     1773         ztmp(:,:) = zmask(:,:) 
     1774         DO ib_bdy = 1, nb_bdy 
     1775            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagv defined only for rims 0 and 1 
     1776               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     1777               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     1778               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10. 
     1779               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 
     1780            END DO 
     1781         END DO 
     1782         CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 
     1783         ztmp(:,:) = zmask(:,:) 
     1784         DO ib_bdy = 1, nb_bdy 
     1785            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! ntreat defined only for rims 0 and 1 
     1786               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     1787               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     1788               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10. 
     1789               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 
     1790            END DO 
     1791         END DO 
     1792         CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 
     1793      END DO 
     1794      CALL iom_close( inum ) 
     1795 
     1796   END SUBROUTINE bdy_meshwri 
     1797    
    17361798   !!================================================================================= 
    17371799END MODULE bdyini 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdylib.F90

    r11191 r11258  
    149149      REAL(wp) ::   zdy_1, zdy_2, zsign_ups 
    150150      REAL(wp), PARAMETER :: zepsilon = 1.e-30                 ! local small value 
    151       REAL(wp), POINTER, DIMENSION(:,:)          :: pmask      ! land/sea mask for field 
    152       REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_xdif ! land/sea mask for x-derivatives 
    153       REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_ydif ! land/sea mask for y-derivatives 
     151      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask      ! land/sea mask for field 
     152      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask_xdif ! land/sea mask for x-derivatives 
     153      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask_ydif ! land/sea mask for y-derivatives 
    154154      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_xdif    ! scale factors for x-derivatives 
    155155      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_ydif    ! scale factors for y-derivatives 
     
    162162      SELECT CASE(igrd) 
    163163         CASE(1) 
    164             pmask      => tmask(:,:,1) 
    165             pmask_xdif => umask(:,:,1) 
    166             pmask_ydif => vmask(:,:,1) 
     164            zmask      => tmask(:,:,1) 
     165            zmask_xdif => umask(:,:,1) 
     166            zmask_ydif => vmask(:,:,1) 
    167167            pe_xdif    => e1u(:,:) 
    168168            pe_ydif    => e2v(:,:) 
     
    170170            ij_offset = 0 
    171171         CASE(2) 
    172             pmask      => umask(:,:,1) 
    173             pmask_xdif => tmask(:,:,1) 
    174             pmask_ydif => fmask(:,:,1) 
     172            zmask      => umask(:,:,1) 
     173            zmask_xdif => tmask(:,:,1) 
     174            zmask_ydif => fmask(:,:,1) 
    175175            pe_xdif    => e1t(:,:) 
    176176            pe_ydif    => e2f(:,:) 
     
    178178            ij_offset = 0 
    179179         CASE(3) 
    180             pmask      => vmask(:,:,1) 
    181             pmask_xdif => fmask(:,:,1) 
    182             pmask_ydif => tmask(:,:,1) 
     180            zmask      => vmask(:,:,1) 
     181            zmask_xdif => fmask(:,:,1) 
     182            zmask_ydif => tmask(:,:,1) 
    183183            pe_xdif    => e1f(:,:) 
    184184            pe_ydif    => e2t(:,:) 
     
    214214         ! 
    215215         ! Calculate scale factors for calculation of spatial derivatives.         
    216          zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1+ii_offset,ijbm1          )         & 
    217         &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1          ,ijbm1+ij_offset) )  
    218          zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2+ii_offset,ijbm2          )         & 
    219         &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2          ,ijbm2+ij_offset) )  
    220          zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )  &  
     216         zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1   +ii_offset,ijbm1             )   & 
     217        &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1             ,ijbm1   +ij_offset) )  
     218         zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2   +ii_offset,ijbm2             )   & 
     219        &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2             ,ijbm2   +ij_offset) )  
     220         zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )   &  
    221221        &      +  (ijbm1-ijbm1jm1) * pe_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) )  
    222          zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1+ii_offset,ijbm1)                  & 
    223         &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1          ,ijbm1+ij_offset) )  
     222         zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1   +ii_offset,ijbm1             )   & 
     223        &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1             ,ijbm1   +ij_offset) )  
    224224         ! make sure scale factors are nonzero 
    225225         if( zey1 .lt. rsmall ) zey1 = zey2 
     
    228228         zey1 = max(zey1,rsmall); zey2 = max(zey2,rsmall);  
    229229         ! 
    230          ! Calculate masks for calculation of spatial derivatives.         
    231          zmask_x = ( abs(iibm1-iibm2) * pmask_xdif(iibm2+ii_offset,ijbm2          )         & 
    232         &          + abs(ijbm1-ijbm2) * pmask_ydif(iibm2          ,ijbm2+ij_offset) )  
    233          zmask_y1 = ( (iibm1-iibm1jm1) * pmask_xdif(iibm1jm1+ii_offset,ijbm1jm1          )  &  
    234         &          +  (ijbm1-ijbm1jm1) * pmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) )  
    235          zmask_y2 = ( (iibm1jp1-iibm1) * pmask_xdif(iibm1+ii_offset,ijbm1)                  & 
    236         &          +  (ijbm1jp1-ijbm1) * pmask_ydif(iibm1          ,ijbm1+ij_offset) )  
     230         ! Calculate masks for calculation of spatial derivatives. 
     231         zmask_x  = ( abs(iibm1-iibm2) * zmask_xdif(iibm2   +ii_offset,ijbm2               )   & 
     232        &           + abs(ijbm1-ijbm2) * zmask_ydif(iibm2             ,ijbm2   +ij_offset) )  
     233         zmask_y1 = ( (iibm1-iibm1jm1) * zmask_xdif(iibm1jm1+ii_offset,ijbm1jm1            )   &  
     234        &          +  (ijbm1-ijbm1jm1) * zmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) )  
     235         zmask_y2 = ( (iibm1jp1-iibm1) * zmask_xdif(iibm1   +ii_offset,ijbm1               )   & 
     236        &          +  (ijbm1jp1-ijbm1) * zmask_ydif(iibm1             ,ijbm1   +ij_offset) )  
    237237 
    238238         ! Calculation of terms required for both versions of the scheme.  
     
    242242         ! Note no rdt factor in expression for zdt because it cancels in the expressions for  
    243243         ! zrx and zry. 
    244          zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1) 
    245          zdx = ( ( phia(iibm1,ijbm1) - phia(iibm2,ijbm2) ) / zex2 ) * zmask_x  
     244         zdt   =     phia(iibm1   ,ijbm1   ) - phib(iibm1   ,ijbm1   ) 
     245         zdx   = ( ( phia(iibm1   ,ijbm1   ) - phia(iibm2   ,ijbm2   ) ) / zex2 ) * zmask_x  
    246246         zdy_1 = ( ( phib(iibm1   ,ijbm1   ) - phib(iibm1jm1,ijbm1jm1) ) / zey1 ) * zmask_y1     
    247          zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1)    ) / zey2 ) * zmask_y2  
     247         zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1   ) ) / zey2 ) * zmask_y2  
    248248         zdy_centred = 0.5 * ( zdy_1 + zdy_2 ) 
    249249!!$         zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1) 
     
    276276           &                    + zwgt * ( phi_ext(jb) - phib(ii,ij) ) ) / ( 1. + zrx )  
    277277         end if 
    278          phia(ii,ij) = phia(ii,ij) * pmask(ii,ij) 
     278         phia(ii,ij) = phia(ii,ij) * zmask(ii,ij) 
    279279      END DO 
    280280      ! 
     
    314314      REAL(wp) ::   zdy_1, zdy_2,  zsign_ups 
    315315      REAL(wp), PARAMETER :: zepsilon = 1.e-30                 ! local small value 
    316       REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field 
    317       REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_xdif ! land/sea mask for x-derivatives 
    318       REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_ydif ! land/sea mask for y-derivatives 
     316      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask      ! land/sea mask for field 
     317      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask_xdif ! land/sea mask for x-derivatives 
     318      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask_ydif ! land/sea mask for y-derivatives 
    319319      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_xdif    ! scale factors for x-derivatives 
    320320      REAL(wp), POINTER, DIMENSION(:,:)          :: pe_ydif    ! scale factors for y-derivatives 
     
    327327      SELECT CASE(igrd) 
    328328         CASE(1) 
    329             pmask      => tmask(:,:,:) 
    330             pmask_xdif => umask(:,:,:) 
    331             pmask_ydif => vmask(:,:,:) 
     329            zmask      => tmask(:,:,:) 
     330            zmask_xdif => umask(:,:,:) 
     331            zmask_ydif => vmask(:,:,:) 
    332332            pe_xdif    => e1u(:,:) 
    333333            pe_ydif    => e2v(:,:) 
     
    335335            ij_offset = 0 
    336336         CASE(2) 
    337             pmask      => umask(:,:,:) 
    338             pmask_xdif => tmask(:,:,:) 
    339             pmask_ydif => fmask(:,:,:) 
     337            zmask      => umask(:,:,:) 
     338            zmask_xdif => tmask(:,:,:) 
     339            zmask_ydif => fmask(:,:,:) 
    340340            pe_xdif    => e1t(:,:) 
    341341            pe_ydif    => e2f(:,:) 
     
    343343            ij_offset = 0 
    344344         CASE(3) 
    345             pmask      => vmask(:,:,:) 
    346             pmask_xdif => fmask(:,:,:) 
    347             pmask_ydif => tmask(:,:,:) 
     345            zmask      => vmask(:,:,:) 
     346            zmask_xdif => fmask(:,:,:) 
     347            zmask_ydif => tmask(:,:,:) 
    348348            pe_xdif    => e1f(:,:) 
    349349            pe_ydif    => e2t(:,:) 
     
    381381            ! 
    382382            ! Calculate scale factors for calculation of spatial derivatives.         
    383             zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1+ii_offset,ijbm1          )         & 
    384            &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1          ,ijbm1+ij_offset) )  
    385             zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2+ii_offset,ijbm2          )         & 
    386            &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2          ,ijbm2+ij_offset) )  
    387             zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )  &  
     383            zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1   +ii_offset,ijbm1             )   & 
     384           &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1             ,ijbm1+ij_offset   ) )  
     385            zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2   +ii_offset,ijbm2             )   & 
     386           &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2             ,ijbm2+ij_offset   ) )  
     387            zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )   &  
    388388           &      +  (ijbm1-ijbm1jm1) * pe_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) )  
    389             zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1+ii_offset,ijbm1)                  & 
    390            &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1          ,ijbm1+ij_offset) )  
     389            zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1   +ii_offset,ijbm1             )   & 
     390           &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1             ,ijbm1+ij_offset   ) )  
    391391            ! make sure scale factors are nonzero 
    392392            if( zey1 .lt. rsmall ) zey1 = zey2 
     
    396396            ! 
    397397            ! Calculate masks for calculation of spatial derivatives.         
    398             zmask_x = ( abs(iibm1-iibm2) * pmask_xdif(iibm2+ii_offset,ijbm2          ,jk)          & 
    399            &          + abs(ijbm1-ijbm2) * pmask_ydif(iibm2          ,ijbm2+ij_offset,jk) )  
    400             zmask_y1 = ( (iibm1-iibm1jm1) * pmask_xdif(iibm1jm1+ii_offset,ijbm1jm1          ,jk)   &  
    401            &          +  (ijbm1-ijbm1jm1) * pmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset,jk) )  
    402             zmask_y2 = ( (iibm1jp1-iibm1) * pmask_xdif(iibm1+ii_offset,ijbm1          ,jk)         & 
    403            &          +  (ijbm1jp1-ijbm1) * pmask_ydif(iibm1          ,ijbm1+ij_offset,jk) )  
     398            zmask_x  = ( abs(iibm1-iibm2) * zmask_xdif(iibm2   +ii_offset,ijbm2             ,jk)   & 
     399           &           + abs(ijbm1-ijbm2) * zmask_ydif(iibm2             ,ijbm2   +ij_offset,jk) )  
     400            zmask_y1 = ( (iibm1-iibm1jm1) * zmask_xdif(iibm1jm1+ii_offset,ijbm1jm1          ,jk)   &  
     401           &          +  (ijbm1-ijbm1jm1) * zmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset,jk) )  
     402            zmask_y2 = ( (iibm1jp1-iibm1) * zmask_xdif(iibm1   +ii_offset,ijbm1             ,jk)   & 
     403           &          +  (ijbm1jp1-ijbm1) * zmask_ydif(iibm1             ,ijbm1   +ij_offset,jk) )  
    404404            ! 
    405405            ! Calculate normal (zrx) and tangential (zry) components of radiation velocities. 
     
    407407            ! Centred derivative is calculated as average of "left" and "right" derivatives for  
    408408            ! this reason.  
    409             zdt = phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk) 
    410             zdx = ( ( phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk) ) / zex2 ) * zmask_x                   
     409            zdt   =     phia(iibm1   ,ijbm1   ,jk) - phib(iibm1   ,ijbm1   ,jk) 
     410            zdx   = ( ( phia(iibm1   ,ijbm1   ,jk) - phia(iibm2   ,ijbm2   ,jk) ) / zex2 ) * zmask_x                   
    411411            zdy_1 = ( ( phib(iibm1   ,ijbm1   ,jk) - phib(iibm1jm1,ijbm1jm1,jk) ) / zey1 ) * zmask_y1   
    412412            zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1   ,ijbm1   ,jk) ) / zey2 ) * zmask_y2       
     
    442442              &                       + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) ) ) / ( 1. + zrx )  
    443443            end if 
    444             phia(ii,ij,jk) = phia(ii,ij,jk) * pmask(ii,ij,jk) 
     444            phia(ii,ij,jk) = phia(ii,ij,jk) * zmask(ii,ij,jk) 
    445445         END DO 
    446446         ! 
     
    469469      !!  
    470470      REAL(wp) ::   zweight 
    471       REAL(wp), POINTER, DIMENSION(:,:,:)      :: pmask         ! land/sea mask for field 
     471      REAL(wp), POINTER, DIMENSION(:,:,:)      :: zmask         ! land/sea mask for field 
    472472      INTEGER  ::   ib, ik   ! dummy loop indices 
    473473      INTEGER  ::   ii, ij   ! 2D addresses 
     
    480480      ! 
    481481      SELECT CASE(igrd) 
    482          CASE(1)   ;   pmask => tmask(:,:,:) 
    483          CASE(2)   ;   pmask => umask(:,:,:) 
    484          CASE(3)   ;   pmask => vmask(:,:,:) 
     482         CASE(1)   ;   zmask => tmask(:,:,:) 
     483         CASE(2)   ;   zmask => umask(:,:,:) 
     484         CASE(3)   ;   zmask => vmask(:,:,:) 
    485485         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' ) 
    486486      END SELECT 
     
    502502            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE 
    503503            DO ik = 1, ipkm1 
    504                IF( pmask(ii1,ij1,ik) /= 0. )   phia(ii,ij,ik) = phia(ii1,ij1,ik)   
     504               IF( zmask(ii1,ij1,ik) /= 0. )   phia(ii,ij,ik) = phia(ii1,ij1,ik)   
    505505            END DO 
    506506         CASE( 9:12 ) 
     
    508508            IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE 
    509509            DO ik = 1, ipkm1 
    510                zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik) 
     510               zweight = zmask(ii1,ij1,ik) + zmask(ii2,ij2,ik) 
    511511               IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) ) / zweight 
    512512            END DO 
     
    516516            IF( ii3 < 1 .OR. ii3 > jpi .OR. ij3 < 1 .OR. ij3 > jpj )   CYCLE 
    517517            DO ik = 1, ipkm1 
    518                zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik) + pmask(ii3,ij3,ik) 
     518               zweight = zmask(ii1,ij1,ik) + zmask(ii2,ij2,ik) + zmask(ii3,ij3,ik) 
    519519               IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) + phia(ii3,ij3,ik) ) / zweight 
    520520            END DO 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DOM/domain.F90

    r10425 r11258  
    101101         CASE( 0 )   ;   WRITE(numout,*) '         (i.e. closed)' 
    102102         CASE( 1 )   ;   WRITE(numout,*) '         (i.e. cyclic east-west)' 
    103          CASE( 2 )   ;   WRITE(numout,*) '         (i.e. equatorial symmetric)' 
     103         CASE( 2 )   ;   WRITE(numout,*) '         (i.e. cyclic north-south)' 
    104104         CASE( 3 )   ;   WRITE(numout,*) '         (i.e. north fold with T-point pivot)' 
    105105         CASE( 4 )   ;   WRITE(numout,*) '         (i.e. cyclic east-west and north fold with T-point pivot)' 
     
    527527      INTEGER ::   inum, ii   ! local integer 
    528528      REAL(wp) ::   zorca_res                     ! local scalars 
    529       REAL(wp) ::   ziglo, zjglo, zkglo, zperio   !   -      - 
     529      REAL(wp) ::   zperio                        !   -      - 
     530      INTEGER, DIMENSION(4) ::   idvar, idimsz    ! size   of dimensions 
    530531      !!---------------------------------------------------------------------- 
    531532      ! 
     
    559560         ! 
    560561      ENDIF 
    561       ! 
    562       CALL iom_get( inum, 'jpiglo', ziglo  )   ;   kpi = NINT( ziglo ) 
    563       CALL iom_get( inum, 'jpjglo', zjglo  )   ;   kpj = NINT( zjglo ) 
    564       CALL iom_get( inum, 'jpkglo', zkglo  )   ;   kpk = NINT( zkglo ) 
     562       ! 
     563      idvar = iom_varid( inum, 'e3t_0', kdimsz = idimsz )   ! use e3t_0, that must exist, to get jp(ijk)glo 
     564      kpi = idimsz(1) 
     565      kpj = idimsz(2) 
     566      kpk = idimsz(3) 
    565567      CALL iom_get( inum, 'jperio', zperio )   ;   kperio = NINT( zperio ) 
    566568      CALL iom_close( inum ) 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynkeg.F90

    r11195 r11258  
    7474      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme  
    7575      ! 
    76       INTEGER  ::   ji, jj, jk, jb           ! dummy loop indices 
    77       INTEGER  ::   igrd, ib_bdy             ! local integers 
     76      INTEGER  ::   ji, jj, jk             ! dummy loop indices 
    7877      REAL(wp) ::   zu, zv                   ! local scalars 
    7978      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke 
    8079      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv  
    81       REAL(wp)  :: zweightu, zweightv 
    82       LOGICAL, DIMENSION(4) :: llsend1, llrecv1  ! indicate how bdy communications are to be carried out 
    8380      !!---------------------------------------------------------------------- 
    8481      ! 
     
    113110            END DO 
    114111         END DO 
    115          ! 
    116          IF (ln_bdy) THEN 
    117             ! Maria Luneva & Fred Wobus: July-2016 
    118             ! compensate for lack of turbulent kinetic energy on liquid bdy points 
    119             DO ib_bdy = 1, nb_bdy 
    120                IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
    121                   igrd = 1           ! compensating null velocity on the bdy 
    122                   DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    123                      ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 1 to jpi 
    124                      jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 1 to jpj 
    125                      IF( ji == 1 .OR. jj == 1 )  CYCLE 
    126                      DO jk = 1, jpkm1 
    127                         zhke(ji,jj,jk) = 0._wp 
    128                         zweightu = umask(ji-1,jj  ,jk) + umask(ji,jj,jk) 
    129                         zweightv = vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk) 
    130                         zu = un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)  +  un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) 
    131                         zv = vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  +  vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) 
    132                         IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zu / (2._wp * zweightu)  
    133                         IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zv / (2._wp * zweightv)  
    134                      END DO 
    135                   END DO 
    136                END IF 
    137             END DO 
    138             ! send jpi-1, jpj-1 and receive 1 used in the computation of the speed tendencies 
    139             llsend1(:) = .false. 
    140             llrecv1(:) = .false. 
    141             DO ib_bdy = 1, nb_bdy 
    142                llsend1(2) = llsend1(2) .OR. ANY(lsend_bdy(ib_bdy,igrd,2,:))   ! send east 
    143                llsend1(4) = llsend1(4) .OR. ANY(lsend_bdy(ib_bdy,igrd,4,:))   ! send north 
    144                llrecv1(1) = llrecv1(1) .OR. ANY(lrecv_bdy(ib_bdy,igrd,1,:))   ! receive west  
    145                llrecv1(3) = llrecv1(3) .OR. ANY(lrecv_bdy(ib_bdy,igrd,3,:))   ! receive south 
    146             END DO 
    147     
    148             IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction 
    149                CALL lbc_lnk( 'bdydyn2d', zhke, 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
    150             END IF 
    151          END IF 
    152          ! 
    153112      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --! 
    154113         DO jk = 1, jpkm1 
     
    168127            END DO 
    169128         END DO 
    170          IF (ln_bdy) THEN 
    171             ! Maria Luneva & Fred Wobus: July-2016 
    172             ! compensate for lack of turbulent kinetic energy on liquid bdy points 
    173             DO ib_bdy = 1, nb_bdy 
    174                IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
    175                   igrd = 1           ! compensation null velocity on land at the bdy 
    176                   DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    177                      ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 1 to jpi 
    178                      jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 1 to jpj 
    179                      IF( ji == 1 .OR. ji == jpi .OR. jj == 1 .OR. jj == jpj )  CYCLE 
    180                      DO jk = 1, jpkm1 
    181                         zhke(ji,jj,jk) = 0._wp 
    182                         zweightu = 8._wp * ( umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) ) & 
    183                              &   + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji  ,jj-1,jk) + umask(ji  ,jj+1,jk) )  
    184                         zweightv = 8._wp * ( vmask(ji  ,jj-1,jk) + vmask(ji  ,jj-1,jk) ) & 
    185                              &   + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj  ,jk) + vmask(ji+1,jj  ,jk) ) 
    186                         zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    & 
    187                            &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  & 
    188                            &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   & 
    189                            &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) 
    190                         zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    & 
    191                            &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  & 
    192                            &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   & 
    193                            &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) 
    194                         IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zu / ( 2._wp * zweightu ) 
    195                         IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zv / ( 2._wp * zweightv ) 
    196                      END DO 
    197                   END DO 
    198                END IF 
    199             END DO 
    200          END IF 
    201129         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. ) 
    202130         ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90

    r11223 r11258  
    465465      ELSE                                      ;   it_offset = 0 
    466466      ENDIF 
    467       IF( PRESENT(kt_offset) )   it_offset = kt_offset 
     467      IF( PRESENT(kt_offset) )      it_offset = kt_offset 
    468468      IF( PRESENT(kit) ) THEN   ;   it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) ) 
    469469      ELSE                      ;   it_offset =         it_offset   * NINT(       rdt            ) 
     
    486486            !       forcing record :    1  
    487487            !                             
    488             ztmp = REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 & 
    489            &       + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday ) 
     488            ztmp =  REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 & 
     489               &  + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday ) 
    490490            sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/)) 
    491491            ! swap at the middle of the year 
     
    517517            !       forcing record :  nmonth  
    518518            !                             
    519             ztmp = REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 & 
    520            &       + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) 
     519            ztmp =  REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 & 
     520           &      + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) 
    521521            imth = nmonth + INT( ztmp ) - COUNT((/llbefore/)) 
    522522            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/)) 
Note: See TracChangeset for help on using the changeset viewer.