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 7339 – NEMO

Changeset 7339


Ignore:
Timestamp:
2016-11-25T16:40:32+01:00 (6 years ago)
Author:
acc
Message:

Ticket #1802. Merged in trunk changes from 6381:6393 and then changes from branches/2016/dev_r6393_NOC_WAD

Location:
branches/2016/dev_NOC_2016/NEMOGCM
Files:
3 deleted
49 edited
23 copied

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_NOC_2016/NEMOGCM/ARCH/arch-macport_osx.fcm

    r5656 r7339  
    5454%CPP               cpp-mp-4.8 
    5555%FC                mpif90  
    56 %FCFLAGS             -fdefault-real-8 -O3 -funroll-all-loops -fcray-pointer  
     56%FCFLAGS             -fdefault-real-8 -O3 -funroll-all-loops -fcray-pointer -ffree-line-length-none 
    5757%FFLAGS              %FCFLAGS 
    5858%LD                  %FC 
  • branches/2016/dev_NOC_2016/NEMOGCM/CONFIG/SHARED/field_def.xml

    r6351 r7339  
    380380      <field_group id="grid_U"   grid_ref="grid_U_2D"> 
    381381         <field id="e3u"          long_name="U-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_U_3D" /> 
     382         <field id="e3u_0"        long_name="Initial U-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_U_3D"/> 
    382383         <field id="utau"         long_name="Wind Stress along i-axis"                               standard_name="surface_downward_x_stress"   unit="N/m2"                            /> 
    383384         <field id="uoce"         long_name="ocean current along i-axis"                             standard_name="sea_water_x_velocity"        unit="m/s"        grid_ref="grid_U_3D" /> 
     
    421422      <field_group id="grid_V"   grid_ref="grid_V_2D"> 
    422423         <field id="e3v"          long_name="V-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_V_3D" /> 
     424         <field id="e3v_0"        long_name="Initial V-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_V_3D"/> 
    423425         <field id="vtau"         long_name="Wind Stress along j-axis"                               standard_name="surface_downward_y_stress"   unit="N/m2"                            /> 
    424426         <field id="voce"         long_name="ocean current along j-axis"                             standard_name="sea_water_y_velocity"        unit="m/s"        grid_ref="grid_V_3D" /> 
  • branches/2016/dev_NOC_2016/NEMOGCM/CONFIG/cfg.txt

    r6140 r7339  
    1111ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 
    1212GYRE OPA_SRC 
     13WAD_TEST_CASES OPA_SRC 
  • branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r6351 r7339  
    145145      CALL iom_put("e3v_0", e3t_0(:,:,:) ) 
    146146      ! 
    147       CALL iom_put( "e3t" , fse3t_n(:,:,:) ) 
    148       CALL iom_put( "e3u" , fse3u_n(:,:,:) ) 
    149       CALL iom_put( "e3v" , fse3v_n(:,:,:) ) 
    150       CALL iom_put( "e3w" , fse3w_n(:,:,:) ) 
     147      CALL iom_put( "e3t" , e3t_n(:,:,:) ) 
     148      CALL iom_put( "e3u" , e3u_n(:,:,:) ) 
     149      CALL iom_put( "e3v" , e3v_n(:,:,:) ) 
     150      CALL iom_put( "e3w" , e3w_n(:,:,:) ) 
    151151      IF( iom_use("e3tdef") )   & 
    152          CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
     152         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    153153 
    154154      CALL iom_put( "ssh" , sshn )                 ! sea surface height 
  • branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r6140 r7339  
    105105      IF( ln_linssh ) THEN          ! Fix in time : set to the reference one for all 
    106106         !       before        !          now          !       after         ! 
    107          ;  gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points 
    108          ;  gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          ! 
    109          ;                     ;   gde3w_n = gde3w_0   !        ---          ! 
     107            gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points 
     108            gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          ! 
     109                                   gde3w_n = gde3w_0   !        ---          ! 
    110110         !                                                                   
    111          ;    e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors 
    112          ;    e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    ! 
    113          ;    e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    ! 
    114          ;                     ;     e3f_n =   e3f_0   !        ---          ! 
    115          ;    e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          ! 
    116          ;   e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          ! 
    117          ;   e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          ! 
     111              e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors 
     112              e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    ! 
     113              e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    ! 
     114                                     e3f_n =   e3f_0   !        ---          ! 
     115              e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          ! 
     116             e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          ! 
     117             e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          ! 
    118118         ! 
    119119         CALL wrk_alloc( jpi,jpj,   z1_hu_0, z1_hv_0 ) 
     
    123123         ! 
    124124         !        before       !          now          !       after         ! 
    125          ;                     ;      ht_n =    ht_0   !                     ! water column thickness 
    126          ;     hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   !  
    127          ;     hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   ! 
    128          ;  r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness 
    129          ;  r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   ! 
     125                                      ht_n =    ht_0   !                     ! water column thickness 
     126               hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   !  
     127               hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   ! 
     128            r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness 
     129            r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   ! 
    130130         ! 
    131131         CALL wrk_dealloc( jpi,jpj,   z1_hu_0, z1_hv_0 ) 
  • branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r6351 r7339  
    874874            ! 
    875875         ELSE                                   !* Initialize at "rest" 
    876             e3t_b(:,:,:) = e3t_0(:,:,:) 
    877             e3t_n(:,:,:) = e3t_0(:,:,:) 
    878             sshn(:,:) = 0.0_wp 
    879  
    880             IF( ln_wd ) THEN 
     876            ! 
     877            IF( ln_wd .AND. ( cp_cfg == 'wad' ) ) THEN 
     878               ! 
     879               CALL wad_istate                  ! WAD test configuration : start from  
     880                                                ! uniform T-S fields and initial ssh slope 
     881               ! needs to be called here and in istate which is called later. 
     882               ! Adjust vertical metrics 
     883               DO jk = 1, jpk 
     884                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 
     885                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
     886                    &            + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk)) 
     887               END DO 
     888               e3t_b(:,:,:) = e3t_n(:,:,:) 
     889               ! 
     890            ELSEIF( ln_wd ) THEN 
     891               ! 
    881892              DO jj = 1, jpj 
    882893                DO ji = 1, jpi 
    883894                  IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN 
    884                      e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1  
    885                      e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1  
    886                      e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1  
     895                     e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 
     896                     e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 
     897                     e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 
    887898                     sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
    888899                     sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     
    891902                ENDDO 
    892903              ENDDO 
     904               ! 
     905            ELSE 
     906               ! 
     907               e3t_b(:,:,:) = e3t_0(:,:,:) 
     908               e3t_n(:,:,:) = e3t_0(:,:,:) 
     909               sshn(:,:) = 0.0_wp 
     910               ! 
    893911            END IF 
    894912 
  • branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r6152 r7339  
    416416               IF(lwp) WRITE(numout,*) '         Depth = rn_bathy read in namelist' 
    417417               zdta(:,:) = rn_bathy 
     418! 
     419               IF( cp_cfg == 'wad' ) THEN 
     420                  SELECT CASE ( jp_cfg ) 
     421                     !                                        ! ==================== 
     422                     CASE ( 1 )                               ! WAD 1 configuration 
     423                        !                                     ! ==================== 
     424                        ! 
     425                        IF(lwp) WRITE(numout,*) 
     426                        IF(lwp) WRITE(numout,*) 'zgr_bat : Closed box with EW linear bottom slope' 
     427                        IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     428                        ! 
     429                        zdta = 1.5_wp 
     430                        DO ji = 10, jpidta 
     431                          zi = MIN(FLOAT(ji - 10)/FLOAT(jpidta - 10), 1.0 ) 
     432                          zdta(ji,:) = MAX(rn_bathy*zi, 1.5)  
     433                          IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 
     434                        END DO 
     435                        !!DO ji = 1, jpidta  
     436                        !!  zi = 1.0-EXP(-0.045*(ji-25.0)**2) 
     437                        !!  zdta(ji,:) = MAX(rn_bathy*zi, 1.5) 
     438                        !!  IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 
     439                        !!END DO 
     440                        zdta(1:2,:) = -2._wp 
     441                        zdta(jpidta-1:jpidta,:) = -2._wp 
     442                        zdta(:,1) = -2._wp 
     443                        zdta(:,jpjdta) = -2._wp 
     444                        zdta(:,1:3) = -2._wp 
     445                        zdta(:,jpjdta-2:jpjdta) = -2._wp 
     446                        !                                     ! ==================== 
     447                     CASE ( 2, 3 )                            ! WAD 2 or 3  configuration 
     448                        !                                     ! ==================== 
     449                        ! 
     450                        IF(lwp) WRITE(numout,*) 
     451                        IF(lwp) WRITE(numout,*) 'zgr_bat : Parobolic EW channel' 
     452                        IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     453                        ! 
     454                        DO ji = 1, jpidta 
     455                          zi = MAX(1.0-FLOAT((ji-25)**2)/484.0, 0.0 ) 
     456                          zi = MAX(1.0-FLOAT((ji-25)**2)/484.0, -2.0 ) 
     457                          zdta(ji,:) = MAX(rn_bathy*zi, -20.0) 
     458                          IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 
     459                        END DO 
     460                        zdta(1:2,:) = -2._wp 
     461                        zdta(jpidta-1:jpidta,:) = -2._wp 
     462                        zdta(:,1) = -2._wp 
     463                        zdta(:,jpjdta) = -2._wp 
     464                        zdta(:,1:3) = -2._wp 
     465                        zdta(:,jpjdta-2:jpjdta) = -2._wp 
     466                        !                                     ! ==================== 
     467                     CASE ( 4 )                               ! WAD 4 configuration 
     468                        !                                     ! ==================== 
     469                        ! 
     470                        IF(lwp) WRITE(numout,*) 
     471                        IF(lwp) WRITE(numout,*) 'zgr_bat : Parobolic bowl' 
     472                        IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     473                        ! 
     474                        DO ji = 1, jpidta 
     475                          zi = MAX(1.0-FLOAT((ji-25)**2)/484.0, -2.0 ) 
     476                        DO jj = 1, jpjdta 
     477                          zj = MAX(1.0-FLOAT((jj-17)**2)/196.0, -2.0 ) 
     478                          zdta(ji,jj) = MAX(rn_bathy*zi*zj, -2.0) 
     479                        END DO 
     480                          IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 
     481                        END DO 
     482                        zdta(1:2,:) = -2._wp 
     483                        zdta(jpidta-1:jpidta,:) = -2._wp 
     484                        zdta(:,1) = -2._wp 
     485                        zdta(:,jpjdta) = -2._wp 
     486                        zdta(:,1:3) = -2._wp 
     487                        zdta(:,jpjdta-2:jpjdta) = -2._wp 
     488                        !                                    ! =========================== 
     489                     CASE ( 5 )                              ! WAD 5 configuration 
     490                        !                                    ! ==================== 
     491                        ! 
     492                        IF(lwp) WRITE(numout,*) 
     493                        IF(lwp) WRITE(numout,*) 'zgr_bat : Double slope with shelf' 
     494                        IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     495                        ! 
     496                        DO ji = 1, jpidta 
     497                          zi = MIN(FLOAT(ji)/FLOAT(jpidta - 5), 1.0 ) 
     498                          zdta(ji,:) = MAX(rn_bathy*zi, 0.5) 
     499                          IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 
     500                        END DO 
     501                        DO ji = jpidta,46,-1 
     502                          zdta(ji,:) = 10.0 
     503                          IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 
     504                        END DO 
     505                        DO ji = 46,20,-1 
     506                          zi = 7.5/25. 
     507                          zdta(ji,:) = MAX(10. - zi*(47.-ji),2.5) 
     508                          IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 
     509                        END DO 
     510                        DO ji = 19,15,-1 
     511                          zdta(ji,:) = 2.5 
     512                          IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 
     513                        END DO 
     514                        DO ji = 15,4,-1 
     515                          zi = 2.0/11.0 
     516                          zdta(ji,:) = MAX(2.5 - zi*(16-ji), 0.5) 
     517                          IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 
     518                        END DO 
     519                        DO ji = 4,1,-1 
     520                          zdta(ji,:) = 0.5 
     521                          IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 
     522                        END DO 
     523                        !                                    ! =========================== 
     524                        zdta(1:2,:) = -4._wp 
     525                        zdta(jpidta-1:jpidta,:) = -4._wp 
     526                        zdta(:,1) = -4._wp 
     527                        zdta(:,jpjdta) = -4._wp 
     528                        zdta(:,1:3) = -4._wp 
     529                        zdta(:,jpjdta-2:jpjdta) = -4._wp 
     530                        !                                    ! =========================== 
     531                     CASE ( 6 )                              ! WAD 6 configuration 
     532                        !                                    ! ==================== 
     533                        ! 
     534                        IF(lwp) WRITE(numout,*) 
     535                        IF(lwp) WRITE(numout,*) 'zgr_bat : Parabolic channel with gaussian ridge' 
     536                        IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     537                        ! 
     538                        DO ji = 1, jpidta 
     539                          zi = MAX(1.0-FLOAT((ji-25)**2)/484.0, -2.0 ) 
     540                          zj = 0.95*MAX(EXP(-1.0*FLOAT((ji-25)**2)/32.0) , 0.0 ) 
     541                          zdta(ji,:) = MAX(rn_bathy*(zi-zj), -2.0) 
     542                          IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 
     543                        END DO 
     544                        zdta(1:2,:) = -4._wp 
     545                        zdta(jpidta-1:jpidta,:) = -4._wp 
     546                        zdta(:,1) = -4._wp 
     547                        zdta(:,jpjdta) = -4._wp 
     548                        zdta(:,1:3) = -4._wp 
     549                        zdta(:,jpjdta-2:jpjdta) = -4._wp 
     550                        !                                    ! =========================== 
     551                     CASE DEFAULT 
     552                        !                                    ! =========================== 
     553                        WRITE(ctmp1,*) 'WAD test with a ', jp_cfg,' option is not coded' 
     554                        ! 
     555                        CALL ctl_stop( ctmp1 ) 
     556                        ! 
     557                  END SELECT 
     558               END IF 
     559               ! 
    418560               IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk 
    419561                  idta(:,:) = jpkm1 
     
    21852327      CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    21862328      ! 
    2187       IF( .NOT.ln_wd ) THEN 
    2188         WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp 
    2189         WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp 
    2190         WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp 
    2191         WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp 
    2192         WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp 
    2193         WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp 
    2194         WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp 
    2195       END IF 
     2329      WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp 
     2330      WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp 
     2331      WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp 
     2332      WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp 
     2333      WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp 
     2334      WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp 
     2335      WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp 
    21962336 
    21972337#if defined key_agrif 
     
    22952435               DO jk = 1, mbathy(ji,jj) 
    22962436                 ! check coordinate is monotonically increasing 
    2297                  IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN 
     2437                 IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN 
    22982438                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    22992439                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    2300                     WRITE(numout,*) 'e3w',e3w_n(ji,jj,:) 
    2301                     WRITE(numout,*) 'e3t',e3t_n(ji,jj,:) 
     2440                    WRITE(numout,*) 'e3w',e3w_0(ji,jj,:) 
     2441                    WRITE(numout,*) 'e3t',e3t_0(ji,jj,:) 
    23022442                    CALL ctl_stop( ctmp1 ) 
    23032443                 ENDIF 
    23042444                 ! and check it has never gone negative 
    2305                  IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN 
     2445                 IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN 
    23062446                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
    23072447                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
    2308                     WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 
    2309                     WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 
     2448                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) 
     2449                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) 
    23102450                    CALL ctl_stop( ctmp1 ) 
    23112451                 ENDIF 
    23122452                 ! and check it never exceeds the total depth 
    2313                  IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     2453                 IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    23142454                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
    23152455                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
    2316                     WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 
     2456                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) 
    23172457                    CALL ctl_stop( ctmp1 ) 
    23182458                 ENDIF 
     
    23212461               DO jk = 1, mbathy(ji,jj)-1 
    23222462                 ! and check it never exceeds the total depth 
    2323                 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     2463                IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    23242464                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
    23252465                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
    2326                     WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 
     2466                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) 
    23272467                    CALL ctl_stop( ctmp1 ) 
    23282468                 ENDIF 
  • branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r6140 r7339  
    3636   USE domvvl          ! varying vertical mesh 
    3737   USE iscplrst        ! ice sheet coupling 
     38   USE wet_dry         ! wetting and drying (needed for wad_istate) 
    3839   ! 
    3940   USE in_out_manager  ! I/O manager 
     
    105106         ELSEIF( cp_cfg == 'gyre' ) THEN          
    106107            CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields 
     108         ELSEIF( cp_cfg == 'wad' ) THEN          
     109            CALL wad_istate                      ! WAD test configuration : start from pre-defined T-S fields and initial ssh slope 
    107110         ELSE                                    ! Initial T-S, U-V fields read in files 
    108111            IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000 
  • branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r6152 r7339  
    432432      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices 
    433433      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars 
    434       LOGICAL  ::   ll_tmp1, ll_tmp2, ll_tmp3            ! local logical variables 
     434      LOGICAL  ::   ll_tmp1, ll_tmp2                     ! local logical variables 
    435435      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
    436436      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy !W/D pressure filter 
     
    438438      ! 
    439439      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
    440       IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
     440      IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    441441      ! 
    442442      IF( kt == nit000 ) THEN 
     
    451451      ENDIF 
    452452      ! 
    453       IF(ln_wd) THEN 
     453      IF( ln_wd ) THEN 
    454454        DO jj = 2, jpjm1 
    455455           DO ji = 2, jpim1  
    456              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))  
    457              ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 
    458              ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
    459                                                        & rn_wdmin1 + rn_wdmin2 
    460  
    461              IF(ll_tmp1.AND.ll_tmp2) THEN 
     456             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
     457                  &    MAX( -bathy(ji,jj)               , -bathy(ji+1,jj) ) .AND.            & 
     458                  &    MAX(   sshn(ji,jj) + bathy(ji,jj),   sshn(ji+1,jj) + bathy(ji+1,jj) ) & 
     459                  &                                                         > rn_wdmin1 + rn_wdmin2 
     460             ll_tmp2 = MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
     461                  &    MAX( -bathy(ji,jj)               , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 
     462 
     463             IF(ll_tmp1) THEN 
    462464               zcpx(ji,jj) = 1.0_wp 
    463                wduflt(ji,jj) = 1.0_wp 
    464              ELSE IF(ll_tmp3) THEN 
    465                ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
    466                zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / & 
    467                            &     (sshn(ji+1,jj) - sshn(ji,jj))) 
    468                wduflt(ji,jj) = 1.0_wp 
     465             ELSE IF(ll_tmp2) THEN 
     466               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     467               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 
     468                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
    469469             ELSE 
    470470               zcpx(ji,jj) = 0._wp 
    471                wduflt(ji,jj) = 0.0_wp 
    472471             END IF 
    473472       
    474              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))  
    475              ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 
    476              ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
    477                                                        & rn_wdmin1 + rn_wdmin2 
    478  
    479              IF(ll_tmp1.AND.ll_tmp2) THEN 
     473             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
     474                  &    MAX( -bathy(ji,jj)               , -bathy(ji,jj+1) ) .AND.            & 
     475                  &    MAX(   sshn(ji,jj) + bathy(ji,jj),   sshn(ji,jj+1) + bathy(ji,jj+1) ) & 
     476                  &                                                         > rn_wdmin1 + rn_wdmin2 
     477             ll_tmp2 = MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
     478                  &    MAX( -bathy(ji,jj)               , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 
     479 
     480             IF(ll_tmp1) THEN 
    480481               zcpy(ji,jj) = 1.0_wp 
    481                wdvflt(ji,jj) = 1.0_wp 
    482              ELSE IF(ll_tmp3) THEN 
    483                ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
    484                zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / & 
    485                            &     (sshn(ji,jj+1) - sshn(ji,jj))) 
    486                wdvflt(ji,jj) = 1.0_wp 
     482             ELSE IF(ll_tmp2) THEN 
     483               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
     484               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 
     485                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) ) 
    487486             ELSE 
    488487               zcpy(ji,jj) = 0._wp 
    489                wdvflt(ji,jj) = 0.0_wp 
    490488             END IF 
    491489           END DO 
    492490        END DO 
    493491        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    494       ENDIF 
    495  
     492      END IF 
    496493 
    497494      ! Surface value 
     
    510507 
    511508 
    512             IF(ln_wd) THEN 
     509            IF( ln_wd ) THEN 
    513510 
    514511              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     
    541538                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 
    542539 
    543                IF(ln_wd) THEN 
     540               IF( ln_wd ) THEN 
    544541                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
    545542                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     
    556553      ! 
    557554      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
    558       IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
     555      IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    559556      ! 
    560557   END SUBROUTINE hpg_sco 
     
    701698      CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    702699      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
    703       IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    704       ! 
    705       ! 
    706       IF(ln_wd) THEN 
     700      IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
     701      ! 
     702      ! 
     703      IF( ln_wd ) THEN 
    707704        DO jj = 2, jpjm1 
    708705           DO ji = 2, jpim1  
    709              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
    710                      & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
    711                      &  > rn_wdmin1 + rn_wdmin2 
    712              ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
    713                                                        & rn_wdmin1 + rn_wdmin2 
     706             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
     707                  &    MAX( -bathy(ji,jj)               , -bathy(ji+1,jj) ) .AND.            & 
     708                  &    MAX(   sshn(ji,jj) + bathy(ji,jj),   sshn(ji+1,jj) + bathy(ji+1,jj) ) & 
     709                  &                                                         > rn_wdmin1 + rn_wdmin2 
     710             ll_tmp2 = MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
     711                  &    MAX( -bathy(ji,jj)               , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 
    714712 
    715713             IF(ll_tmp1) THEN 
    716714               zcpx(ji,jj) = 1.0_wp 
    717715             ELSE IF(ll_tmp2) THEN 
    718                ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
    719                zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
    720                            &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     716               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     717               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 
     718                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
    721719             ELSE 
    722720               zcpx(ji,jj) = 0._wp 
    723721             END IF 
    724722       
    725              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
    726                      & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
    727                      &  > rn_wdmin1 + rn_wdmin2 
    728              ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
    729                                                        & rn_wdmin1 + rn_wdmin2 
     723             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
     724                  &    MAX( -bathy(ji,jj)               , -bathy(ji,jj+1) ) .AND.            & 
     725                  &    MAX(   sshn(ji,jj) + bathy(ji,jj),   sshn(ji,jj+1) + bathy(ji,jj+1) ) & 
     726                  &                                                         > rn_wdmin1 + rn_wdmin2 
     727             ll_tmp2 = MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
     728                  &    MAX( -bathy(ji,jj)               , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 
    730729 
    731730             IF(ll_tmp1) THEN 
    732731               zcpy(ji,jj) = 1.0_wp 
    733732             ELSE IF(ll_tmp2) THEN 
    734                ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
    735                zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
    736                            &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     733               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
     734               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 
     735                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) ) 
    737736             ELSE 
    738737               zcpy(ji,jj) = 0._wp 
     
    741740        END DO 
    742741        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    743       ENDIF 
    744  
     742      END IF 
    745743 
    746744      IF( kt == nit000 ) THEN 
     
    913911            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 
    914912            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 
    915             IF(ln_wd) THEN 
     913            IF( ln_wd ) THEN 
    916914              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
    917915              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     
    936934                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    937935                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj) 
    938                IF(ln_wd) THEN 
     936               IF( ln_wd ) THEN 
    939937                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
    940938                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     
    950948      CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    951949      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
    952       IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
     950      IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    953951      ! 
    954952   END SUBROUTINE hpg_djc 
     
    987985      CALL wrk_alloc( jpi,jpj,jpk,   zdept, zrhh ) 
    988986      CALL wrk_alloc( jpi,jpj,       zsshu_n, zsshv_n ) 
    989       IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
     987      IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    990988      ! 
    991989      IF( kt == nit000 ) THEN 
     
    1000998      IF( ln_linssh )   znad = 0._wp 
    1001999 
    1002       IF(ln_wd) THEN 
     1000      IF( ln_wd ) THEN 
    10031001        DO jj = 2, jpjm1 
    10041002           DO ji = 2, jpim1  
    1005              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
    1006                      & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
    1007                      &  > rn_wdmin1 + rn_wdmin2 
    1008              ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
    1009                                                        & rn_wdmin1 + rn_wdmin2 
     1003             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
     1004                  &    MAX( -bathy(ji,jj)               , -bathy(ji+1,jj) ) .AND.            & 
     1005                  &    MAX(   sshn(ji,jj) + bathy(ji,jj),   sshn(ji+1,jj) + bathy(ji+1,jj) ) & 
     1006                  &                                                         > rn_wdmin1 + rn_wdmin2 
     1007             ll_tmp2 = MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
     1008                  &    MAX( -bathy(ji,jj)               , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 
    10101009 
    10111010             IF(ll_tmp1) THEN 
    10121011               zcpx(ji,jj) = 1.0_wp 
    10131012             ELSE IF(ll_tmp2) THEN 
    1014                ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
    1015                zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
    1016                            &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     1013               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     1014               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 
     1015                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
    10171016             ELSE 
    10181017               zcpx(ji,jj) = 0._wp 
    10191018             END IF 
    10201019       
    1021              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
    1022                      & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
    1023                      &  > rn_wdmin1 + rn_wdmin2 
    1024              ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
    1025                                                        & rn_wdmin1 + rn_wdmin2 
    1026  
    1027              IF(ll_tmp1.OR.ll_tmp2) THEN 
     1020             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
     1021                  &    MAX( -bathy(ji,jj)               , -bathy(ji,jj+1) ) .AND.            & 
     1022                  &    MAX(   sshn(ji,jj) + bathy(ji,jj),   sshn(ji,jj+1) + bathy(ji,jj+1) ) & 
     1023                  &                                                         > rn_wdmin1 + rn_wdmin2 
     1024             ll_tmp2 = MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
     1025                  &    MAX( -bathy(ji,jj)               , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 
     1026 
     1027             IF(ll_tmp1) THEN 
    10281028               zcpy(ji,jj) = 1.0_wp 
    10291029             ELSE IF(ll_tmp2) THEN 
    1030                ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
    1031                zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
    1032                            &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     1030               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
     1031               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 
     1032                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) ) 
    10331033             ELSE 
    10341034               zcpy(ji,jj) = 0._wp 
     
    10371037        END DO 
    10381038        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    1039       ENDIF 
     1039      END IF 
    10401040 
    10411041      ! Clean 3-D work arrays 
     
    12211221                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
    12221222               ENDIF 
    1223                IF(ln_wd) THEN 
     1223               IF( ln_wd ) THEN 
    12241224                  zdpdx1 = zdpdx1 * zcpx(ji,jj) 
    12251225                  zdpdx2 = zdpdx2 * zcpx(ji,jj) 
     
    12801280                  zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
    12811281               ENDIF 
    1282                IF(ln_wd) THEN 
     1282               IF( ln_wd ) THEN 
    12831283                  zdpdy1 = zdpdy1 * zcpy(ji,jj) 
    12841284                  zdpdy2 = zdpdy2 * zcpy(ji,jj) 
     
    12951295      CALL wrk_dealloc( jpi,jpj,jpk,   zdept, zrhh ) 
    12961296      CALL wrk_dealloc( jpi,jpj,       zsshu_n, zsshv_n ) 
    1297       IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
     1297      IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    12981298      ! 
    12991299   END SUBROUTINE hpg_prj 
  • branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r6140 r7339  
    205205      ENDIF 
    206206      !                          ! Control of surface pressure gradient scheme options 
    207       ;                              nspg =  np_NO    ;   ioptio = 0 
     207                                     nspg =  np_NO    ;   ioptio = 0 
    208208      IF( ln_dynspg_exp ) THEN   ;   nspg =  np_EXP   ;   ioptio = ioptio + 1   ;   ENDIF 
    209209      IF( ln_dynspg_ts  ) THEN   ;   nspg =  np_TS    ;   ioptio = ioptio + 1   ;   ENDIF 
  • branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r6152 r7339  
    156156      REAL(wp), POINTER, DIMENSION(:,:) :: zhf 
    157157      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef. 
    158       REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1           ! Wetting/Dying velocity filter coef. 
    159158      !!---------------------------------------------------------------------- 
    160159      ! 
     
    168167      CALL wrk_alloc( jpi,jpj,   zsshu_a, zsshv_a                  ) 
    169168      CALL wrk_alloc( jpi,jpj,   zhf ) 
    170       IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
     169      IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
    171170      ! 
    172171      zmdi=1.e+20                               !  missing data indicator for masking 
     
    374373      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient 
    375374        IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters 
    376           wduflt1(:,:) = 1.0_wp 
    377           wdvflt1(:,:) = 1.0_wp 
    378           DO jj = 2, jpjm1 
    379              DO ji = 2, jpim1 
    380                 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))  & 
    381                         & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj))   & 
    382                         &  > rn_wdmin1 + rn_wdmin2 
    383                 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
    384                         &                                   + rn_wdmin1 + rn_wdmin2 
     375           DO jj = 2, jpjm1 
     376              DO ji = 2, jpim1  
     377                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
     378                     &    MAX( -bathy(ji,jj)               , -bathy(ji+1,jj) ) .AND.            & 
     379                     &    MAX(   sshn(ji,jj) + bathy(ji,jj),   sshn(ji+1,jj) + bathy(ji+1,jj) ) & 
     380                     &                                                         > rn_wdmin1 + rn_wdmin2 
     381                ll_tmp2 = MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
     382                     &    MAX( -bathy(ji,jj)               , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 
     383    
    385384                IF(ll_tmp1) THEN 
    386                   zcpx(ji,jj)    = 1.0_wp 
    387                 ELSEIF(ll_tmp2) THEN 
    388                    ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here 
    389                   zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 
    390                         &          /(sshn(ji+1,jj) - sshn(ji,jj))) 
     385                  zcpx(ji,jj) = 1.0_wp 
     386                ELSE IF(ll_tmp2) THEN 
     387                  ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     388                  zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 
     389                              &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
    391390                ELSE 
    392                   zcpx(ji,jj)    = 0._wp 
    393                   wduflt1(ji,jj) = 0.0_wp 
     391                  zcpx(ji,jj) = 0._wp 
    394392                END IF 
    395  
    396                 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
    397                         & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1))   & 
    398                         &  > rn_wdmin1 + rn_wdmin2 
    399                 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
    400                         &                                   + rn_wdmin1 + rn_wdmin2 
     393          
     394                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
     395                     &    MAX( -bathy(ji,jj)               , -bathy(ji,jj+1) ) .AND.            & 
     396                     &    MAX(   sshn(ji,jj) + bathy(ji,jj),   sshn(ji,jj+1) + bathy(ji,jj+1) ) & 
     397                     &                                                         > rn_wdmin1 + rn_wdmin2 
     398                ll_tmp2 = MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
     399                     &    MAX( -bathy(ji,jj)               , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 
     400    
    401401                IF(ll_tmp1) THEN 
    402                    zcpy(ji,jj)    = 1.0_wp 
    403                 ELSEIF(ll_tmp2) THEN 
    404                    ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here 
    405                   zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 
    406                         &          /(sshn(ji,jj+1) - sshn(ji,jj))) 
     402                  zcpy(ji,jj) = 1.0_wp 
     403                ELSE IF(ll_tmp2) THEN 
     404                  ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
     405                  zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 
     406                              &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) ) 
    407407                ELSE 
    408                   zcpy(ji,jj)    = 0._wp 
    409                   wdvflt1(ji,jj) = 0.0_wp 
    410                 ENDIF 
    411  
    412              END DO 
     408                  zcpy(ji,jj) = 0._wp 
     409                END IF 
     410              END DO 
    413411           END DO 
    414  
    415            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    416  
     412  
    417413           DO jj = 2, jpjm1 
    418414              DO ji = 2, jpim1 
    419                  zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
    420                         &                        * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 
    421                  zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
    422                         &                        * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 
     415                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
     416                        &                        * r1_e1u(ji,jj) * zcpx(ji,jj) 
     417                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
     418                        &                        * r1_e2v(ji,jj) * zcpy(ji,jj) 
    423419              END DO 
    424420           END DO 
     
    567563      ENDIF 
    568564 
    569       IF( ln_wd ) THEN      !preserve the positivity of water depth 
    570                           !ssh[b,n,a] should have already been processed for this 
    571          sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 
    572          sshb_e(:,:)  = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 
    573       ENDIF 
    574565      ! 
    575566      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
     
    646637            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
    647638            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
    648             IF( ln_wd ) THEN 
    649               zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 
    650               zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 
    651             END IF 
    652639         ELSE 
    653640            zhup2_e (:,:) = hu_n(:,:) 
     
    701688            END DO 
    702689         END DO 
     690 
    703691         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:) 
    704          IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))  
     692          
    705693         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
    706694 
     
    749737         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) & 
    750738          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
     739 
    751740         IF( ln_wd ) THEN                   ! Calculating and applying W/D gravity filters 
    752            wduflt1(:,:) = 1._wp 
    753            wdvflt1(:,:) = 1._wp 
    754741           DO jj = 2, jpjm1 
    755               DO ji = 2, jpim1 
    756                  ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 
    757                         & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) )    & 
    758                         &                                  > rn_wdmin1 + rn_wdmin2 
    759                  ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 
    760                         &                                  + rn_wdmin1 + rn_wdmin2 
    761                  IF(ll_tmp1) THEN 
    762                     zcpx(ji,jj) = 1._wp 
    763                  ELSE IF(ll_tmp2) THEN 
    764                     ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen here 
    765                     zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 
    766                         &             / (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)) ) 
    767                  ELSE 
    768                     zcpx(ji,jj)    = 0._wp 
    769                     wduflt1(ji,jj) = 0._wp 
    770                  END IF 
    771  
    772                  ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 
    773                         & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) )    & 
    774                         &                                  > rn_wdmin1 + rn_wdmin2 
    775                  ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 
    776                         &                                  + rn_wdmin1 + rn_wdmin2 
    777                  IF(ll_tmp1) THEN 
    778                     zcpy(ji,jj) = 1._wp 
    779                  ELSE IF(ll_tmp2) THEN 
    780                     ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen here 
    781                     zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 
    782                         &             / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)) ) 
    783                  ELSE 
    784                     zcpy(ji,jj)    = 0._wp 
    785                     wdvflt1(ji,jj) = 0._wp 
    786                  END IF 
     742              DO ji = 2, jpim1  
     743                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                & 
     744                     &    MAX(   -bathy(ji,jj)               ,   -bathy(ji+1,jj) ) .AND.            & 
     745                     &    MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) ) & 
     746                     &                                                             > rn_wdmin1 + rn_wdmin2 
     747                ll_tmp2 = MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                & 
     748                     &    MAX(   -bathy(ji,jj)               ,   -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 
     749    
     750                IF(ll_tmp1) THEN 
     751                  zcpx(ji,jj) = 1.0_wp 
     752                ELSE IF(ll_tmp2) THEN 
     753                  ! no worries about  zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj) = 0, it won't happen ! here 
     754                  zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +    bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 
     755                              &    / (zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj)) ) 
     756                ELSE 
     757                  zcpx(ji,jj) = 0._wp 
     758                END IF 
     759          
     760                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
     761                     &    MAX(   -bathy(ji,jj)               ,   -bathy(ji,jj+1) ) .AND.            & 
     762                     &    MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) ) & 
     763                     &                                                             > rn_wdmin1 + rn_wdmin2 
     764                ll_tmp2 = MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
     765                     &    MAX(   -bathy(ji,jj)               ,   -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 
     766    
     767                IF(ll_tmp1) THEN 
     768                  zcpy(ji,jj) = 1.0_wp 
     769                ELSE IF(ll_tmp2) THEN 
     770                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here 
     771                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +    bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 
     772                              &    / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  )) ) 
     773                ELSE 
     774                  zcpy(ji,jj) = 0._wp 
     775                END IF 
    787776              END DO 
    788             END DO 
    789             CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    790          ENDIF 
     777           END DO 
     778         END IF 
    791779         ! 
    792780         ! Compute associated depths at U and V points: 
     
    806794            END DO 
    807795 
    808             IF( ln_wd ) THEN 
    809               zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 
    810               zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 
    811             END IF 
    812  
    813796         ENDIF 
    814797         ! 
     
    888871                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 
    889872                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 
    890                  zwx(ji,jj) = zu_spg * zcpx(ji,jj)  
    891                  zwy(ji,jj) = zv_spg * zcpy(ji,jj) 
     873                 zwx(ji,jj) = zu_spg * zcpx(ji,jj) * wdmask(ji,jj) * wdmask(ji+1, jj)  
     874                 zwy(ji,jj) = zv_spg * zcpy(ji,jj) * wdmask(ji,jj) * wdmask(ji, jj+1) 
    892875              END DO 
    893876           END DO 
     
    927910               DO ji = fs_2, fs_jpim1   ! vector opt. 
    928911 
    929                   IF( ln_wd ) THEN 
    930                     zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 
    931                     zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 
    932                   ELSE 
    933                     zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
    934                     zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
    935                   END IF 
     912                  zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
     913                  zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
    936914                  zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 
    937915                  zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 
     
    953931         ! 
    954932         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only) 
    955             IF( ln_wd ) THEN 
    956               hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 
    957               hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 
    958             ELSE 
    959               hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
    960               hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
    961             END IF 
     933            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     934            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
    962935            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 
    963936            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 
     
    1024997      ! 
    1025998      ! Update barotropic trend: 
    1026       IF( ln_dynadv_vec .OR. ln_linssh ) THEN 
    1027          DO jk=1,jpkm1 
    1028             ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 
    1029             va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 
    1030          END DO 
     999      IF(ln_wd) THEN 
     1000        IF( ln_dynadv_vec .OR. ln_linssh ) THEN 
     1001           DO jk=1,jpkm1 
     1002              ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b * wdmask(:,:) 
     1003              va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b * wdmask(:,:) 
     1004           END DO 
     1005        ELSE 
     1006           ! At this stage, ssha has been corrected: compute new depths at velocity points 
     1007           DO jj = 1, jpjm1 
     1008              DO ji = 1, jpim1      ! NO Vector Opt. 
     1009                 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj) & 
     1010                    &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)    & 
     1011                    &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 
     1012                 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj) & 
     1013                    &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )    & 
     1014                    &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 
     1015              END DO 
     1016           END DO 
     1017           CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 
     1018           ! 
     1019           DO jk=1,jpkm1 
     1020              ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b * wdmask(:,:) 
     1021              va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b * wdmask(:,:) 
     1022           END DO 
     1023           ! Save barotropic velocities not transport: 
     1024           ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 
     1025           va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
     1026        ENDIF 
    10311027      ELSE 
    1032          ! At this stage, ssha has been corrected: compute new depths at velocity points 
    1033          DO jj = 1, jpjm1 
    1034             DO ji = 1, jpim1      ! NO Vector Opt. 
    1035                zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj) & 
    1036                   &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)    & 
    1037                   &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 
    1038                zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj) & 
    1039                   &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )    & 
    1040                   &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 
    1041             END DO 
    1042          END DO 
    1043          CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 
    1044          ! 
    1045          DO jk=1,jpkm1 
    1046             ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 
    1047             va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 
    1048          END DO 
    1049          ! Save barotropic velocities not transport: 
    1050          ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 
    1051          va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
    1052       ENDIF 
     1028        IF( ln_dynadv_vec .OR. ln_linssh ) THEN 
     1029           DO jk=1,jpkm1 
     1030              ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 
     1031              va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 
     1032           END DO 
     1033        ELSE 
     1034           ! At this stage, ssha has been corrected: compute new depths at velocity points 
     1035           DO jj = 1, jpjm1 
     1036              DO ji = 1, jpim1      ! NO Vector Opt. 
     1037                 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj) & 
     1038                    &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)    & 
     1039                    &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 
     1040                 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj) & 
     1041                    &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )    & 
     1042                    &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 
     1043              END DO 
     1044           END DO 
     1045           CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 
     1046           ! 
     1047           DO jk=1,jpkm1 
     1048              ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 
     1049              va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 
     1050           END DO 
     1051           ! Save barotropic velocities not transport: 
     1052           ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 
     1053           va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
     1054        ENDIF 
     1055 
     1056      END IF 
    10531057      ! 
    10541058      DO jk = 1, jpkm1 
     
    10861090      CALL wrk_dealloc( jpi,jpj,   zsshu_a, zsshv_a                                   ) 
    10871091      CALL wrk_dealloc( jpi,jpj,   zhf ) 
    1088       IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
     1092      IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 
    10891093      ! 
    10901094      IF ( ln_diatmb ) THEN 
  • branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r6152 r7339  
    8888      ENDIF 
    8989      ! 
    90       CALL div_hor( kt )                              ! Horizontal divergence 
    91       ! 
    92       z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
     90      z2dt = 2._wp * rdt                          ! set time step size (Euler/Leapfrog) 
    9391      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
     92      zcoef = 0.5_wp * r1_rau0 
    9493 
    9594      !                                           !------------------------------! 
    9695      !                                           !   After Sea Surface Height   ! 
    9796      !                                           !------------------------------! 
     97      IF(ln_wd) THEN 
     98         CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 
     99      ENDIF 
     100 
     101      CALL div_hor( kt )                               ! Horizontal divergence 
     102      ! 
    98103      zhdiv(:,:) = 0._wp 
    99104      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
     
    104109      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    105110      !  
    106       zcoef = 0.5_wp * r1_rau0 
    107  
    108       IF(ln_wd) CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 
    109  
    110111      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
    111112 
  • branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90

    r6152 r7339  
    3333   !! --------------------------------------------------------------------- 
    3434 
    35    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wduflt, wdvflt !: u- and v- filter 
    3635   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask         !: u- and v- limiter  
    3736 
     
    4645   PUBLIC   wad_lmt                   ! routine called by sshwzv.F90 
    4746   PUBLIC   wad_lmt_bt                ! routine called by dynspg_ts.F90 
     47   PUBLIC   wad_istate                ! routine called by istate.F90 and domvvl.F90 
    4848 
    4949   !! * Substitutions 
     
    8787 
    8888      IF(ln_wd) THEN 
    89          ALLOCATE( wduflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr ) 
     89         ALLOCATE( wdmask(jpi,jpj), STAT=ierr ) 
    9090         IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 
    9191      ENDIF 
     
    145145        ! Horizontal Flux in u and v direction 
    146146        DO jk = 1, jpkm1   
    147            DO jj = 1, jpjm1 
    148               DO ji = 1, jpim1 
     147           DO jj = 1, jpj 
     148              DO ji = 1, jpi 
    149149                 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 
    150150                 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 
     
    156156        zflxv(:,:) = zflxv(:,:) * e1v(:,:) 
    157157        
    158         DO jj = 2, jpjm1 
    159            DO ji = 2, jpim1  
     158        wdmask(:,:) = 1 
     159        DO jj = 2, jpj 
     160           DO ji = 2, jpi  
    160161 
    161162             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells 
     
    168169 
    169170              zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 
    170               IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary 
     171              IF(zdep2 .le. 0._wp) THEN  !add more safty, but not necessary 
    171172                !zdep2 = 0._wp 
    172173                sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     174                wdmask(ji,jj) = 0._wp 
    173175              END IF 
    174176           ENDDO 
     
    183185           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 
    184186           
    185            DO jj = 2, jpjm1 
    186               DO ji = 2, jpim1  
     187           DO jj = 2, jpj 
     188              DO ji = 2, jpi  
    187189         
    188                  wdmask(ji,jj) = 0 
    189190                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE  
    190191                 IF(bathy(ji,jj) > zdepwd) CYCLE 
     
    202203                 IF(zdep1 > zdep2) THEN 
    203204                   zflag = 1 
    204                    wdmask(ji, jj) = 1 
     205                   wdmask(ji, jj) = 0 
    205206                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
    206207                   zcoef = max(zcoef, 0._wp) 
     
    209210                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 
    210211                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef 
    211                    IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef 
     212                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 
    212213                 END IF 
    213214              END DO ! ji loop 
     
    231232        CALL lbc_lnk( un, 'U', -1. ) 
    232233        CALL lbc_lnk( vn, 'V', -1. ) 
     234      ! 
     235        un_b(:,:) = un_b(:,:) * zwdlmtu(:, :) 
     236        vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :) 
     237        CALL lbc_lnk( un_b, 'U', -1. ) 
     238        CALL lbc_lnk( vn_b, 'V', -1. ) 
    233239        
    234240        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 
     
    291297        zflxp(:,:)   = 0._wp 
    292298        zflxn(:,:)   = 0._wp 
    293         !zflxu(:,:)   = 0._wp 
    294         !zflxv(:,:)   = 0._wp 
    295299 
    296300        zwdlmtu(:,:)  = 1._wp 
     
    299303        ! Horizontal Flux in u and v direction 
    300304        
    301         !zflxu(:,:) = zflxu(:,:) * e2u(:,:) 
    302         !zflxv(:,:) = zflxv(:,:) * e1v(:,:) 
    303         
    304         DO jj = 2, jpjm1 
    305            DO ji = 2, jpim1  
     305        DO jj = 2, jpj 
     306           DO ji = 2, jpi  
    306307 
    307308             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells 
     
    314315 
    315316              zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 
    316               IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary 
    317                 !zdep2 = 0._wp 
    318                sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
    319               END IF 
    320317           ENDDO 
    321318        END DO 
     
    329326           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 
    330327           
    331            DO jj = 2, jpjm1 
    332               DO ji = 2, jpim1  
     328           DO jj = 2, jpj 
     329              DO ji = 2, jpi  
    333330         
    334                  wdmask(ji,jj) = 0 
    335331                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE  
    336332                 IF(bathy(ji,jj) > zdepwd) CYCLE 
     
    349345                 IF(zdep1 > zdep2) THEN 
    350346                   zflag = 1 
    351                    !wdmask(ji, jj) = 1 
    352347                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
    353348                   zcoef = max(zcoef, 0._wp) 
     
    356351                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 
    357352                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef 
    358                    IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef 
     353                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 
    359354                 END IF 
    360355              END DO ! ji loop 
     
    379374        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 
    380375        
    381         !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field) 
    382         !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field) 
    383376        ! 
    384377        ! 
     
    390383      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt') 
    391384   END SUBROUTINE wad_lmt_bt 
     385 
     386   SUBROUTINE wad_istate 
     387      !!---------------------------------------------------------------------- 
     388      !!                   ***  ROUTINE wad_istate  *** 
     389      !!  
     390      !! ** Purpose :   Initialization of the dynamics and tracers for WAD test 
     391      !!      configurations (channels or bowls with initial ssh gradients) 
     392      !! 
     393      !! ** Method  : - set temperature field 
     394      !!              - set salinity field 
     395      !!              - set ssh slope (needs to be repeated in domvvl_rst_init to 
     396      !!                               set vertical metrics ) 
     397      !!---------------------------------------------------------------------- 
     398      ! 
     399      INTEGER  ::   ji, jj            ! dummy loop indices 
     400      REAL(wp) ::   zi, zj 
     401      !!---------------------------------------------------------------------- 
     402      ! 
     403      ! Uniform T & S in all test cases 
     404      tsn(:,:,:,jp_tem) = 10._wp 
     405      tsb(:,:,:,jp_tem) = 10._wp 
     406      tsn(:,:,:,jp_sal) = 35._wp 
     407      tsb(:,:,:,jp_sal) = 35._wp 
     408      SELECT CASE ( jp_cfg )  
     409         !                                        ! ==================== 
     410         CASE ( 1 )                               ! WAD 1 configuration 
     411            !                                     ! ==================== 
     412            ! 
     413            IF(lwp) WRITE(numout,*) 
     414            IF(lwp) WRITE(numout,*) 'istate_wad : Closed box with EW linear bottom slope' 
     415            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     416            ! 
     417            do ji = 1,jpi 
     418             sshn(ji,:) = ( -5.5_wp + 5.5_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     419            end do 
     420            !                                     ! ==================== 
     421         CASE ( 2 )                               ! WAD 2 configuration 
     422            !                                     ! ==================== 
     423            ! 
     424            IF(lwp) WRITE(numout,*) 
     425            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, mid-range initial ssh slope' 
     426            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     427            ! 
     428            do ji = 1,jpi 
     429             sshn(ji,:) = ( -5.5_wp + 3.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     430            end do 
     431            !                                     ! ==================== 
     432         CASE ( 3 )                               ! WAD 3 configuration 
     433            !                                     ! ==================== 
     434            ! 
     435            IF(lwp) WRITE(numout,*) 
     436            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, extreme initial ssh slope'  
     437            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     438            ! 
     439            do ji = 1,jpi 
     440             sshn(ji,:) = ( -7.5_wp + 6.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     441            end do 
     442 
     443            ! 
     444            !                                     ! ==================== 
     445         CASE ( 4 )                               ! WAD 4 configuration 
     446            !                                     ! ==================== 
     447            ! 
     448            IF(lwp) WRITE(numout,*) 
     449            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic bowl, mid-range initial ssh slope'  
     450            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     451            ! 
     452            DO ji = 1, jpi 
     453               zi = MAX(1.0-FLOAT((mig(ji)-25)**2)/400.0, 0.0 ) 
     454               DO jj = 1, jpj 
     455                  zj = MAX(1.0-FLOAT((mjg(jj)-17)**2)/144.0, 0.0 ) 
     456                  sshn(ji,jj) = -8.5_wp + 8.5_wp*zi*zj 
     457               END DO 
     458            END DO 
     459 
     460            ! 
     461            !                                    ! =========================== 
     462         CASE ( 5 )                              ! WAD 5 configuration 
     463            !                                    ! ==================== 
     464            ! 
     465            IF(lwp) WRITE(numout,*) 
     466            IF(lwp) WRITE(numout,*) 'istate_wad : Double slope with shelf' 
     467            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     468            ! 
     469            ! Needed rn_wdmin2 increased to 0.01 for this case? 
     470            do ji = 1,jpi 
     471             sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     472            end do 
     473 
     474            ! 
     475            !                                     ! =========================== 
     476         CASE ( 6 )                               ! WAD 6 configuration 
     477            !                                     ! ==================== 
     478            ! 
     479            IF(lwp) WRITE(numout,*) 
     480            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel with gaussian ridge'  
     481            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     482            ! 
     483            do ji = 1,jpi 
     484             !6a 
     485             sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     486             !Some variations in initial slope that have been tested 
     487             !6b 
     488             !sshn(ji,:) = ( -5.5_wp + 6.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     489             !6c 
     490             !sshn(ji,:) = ( -5.5_wp + 7.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     491             !6d 
     492             !sshn(ji,:) = ( -4.5_wp + 8.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     493            end do 
     494 
     495            ! 
     496            !                                    ! =========================== 
     497         CASE DEFAULT                            ! NONE existing configuration 
     498            !                                    ! =========================== 
     499            WRITE(ctmp1,*) 'WAD test with a ', jp_cfg,' option is not coded' 
     500            ! 
     501            CALL ctl_stop( ctmp1 ) 
     502            ! 
     503      END SELECT 
     504      ! 
     505      ! Apply minimum wetdepth criterion 
     506      ! 
     507      do jj = 1,jpj 
     508         do ji = 1,jpi 
     509            IF( bathy(ji,jj) + sshn(ji,jj) < rn_wdmin1 ) THEN 
     510               sshn(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - bathy(ji,jj) ) 
     511            ENDIF 
     512         end do 
     513      end do 
     514      sshb = sshn 
     515      ssha = sshn 
     516      ! 
     517   END SUBROUTINE wad_istate 
     518 
     519   !!===================================================================== 
    392520END MODULE wet_dry 
  • branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90

    r6140 r7339  
    149149            DO jj = 2, jpjm1 
    150150               DO ji = fs_2, fs_jpim1   ! vector opt. 
    151                   ! total intermediate advective trends 
     151                  !                             ! total intermediate advective trends 
    152152                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    153153                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    154                      &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    155                   ! update and guess with monotonic sheme 
    156 !!gm why tmask added in the two following lines ???    the mask is done in tranxt ! 
    157                   pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra   * tmask(ji,jj,jk) 
    158                   zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + p2dt * ztra ) * tmask(ji,jj,jk) 
     154                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     155                  !                             ! update and guess with monotonic sheme 
     156                  pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 
     157                  zwi(ji,jj,jk)    = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
    159158               END DO 
    160159            END DO 
     
    163162         !                 
    164163         IF( l_trd )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
    165             ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:) 
     164            ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:) 
    166165         END IF 
    167166         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     
    364363      ! 
    365364      CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav ) 
    366       CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 
     365      CALL wrk_alloc( jpi,jpj,jpk,         zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 
    367366      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs ) 
    368367      ! 
     
    436435                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    437436                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    438                      &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     437                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   ) * r1_e1e2t(ji,jj) 
    439438                  !                             ! update and guess with monotonic sheme 
    440                   pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra 
    441                   zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + p2dt * ztra ) * tmask(ji,jj,jk) 
     439                  pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 
     440                  zwi(ji,jj,jk)    = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
    442441               END DO 
    443442            END DO 
     
    488487         zwz_sav(:,:,:)   = zwz(:,:,:) 
    489488         ztrs   (:,:,:,1) = ptb(:,:,:,jn) 
     489         ztrs   (:,:,1,2) = ptb(:,:,1,jn) 
     490         ztrs   (:,:,1,3) = ptb(:,:,1,jn) 
    490491         zwzts  (:,:,:)   = 0._wp 
    491492         ! 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/DMP_TOOLS/src/zoom.F90

    r4739 r7339  
    2929      NAMELIST/nam_zoom_dmp/lzoom_n,lzoom_e,lzoom_w,lzoom_s 
    3030      !!---------------------------------------------------------------------- 
    31       ! 
    32       IF( nn_timing == 1 )  CALL timing_start( 'dtacof_zoom') 
    33       ! 
    34        
     31 
    3532      ! Read namelist 
    3633      OPEN( UNIT=numnam, FILE='namelist', FORM='FORMATTED', STATUS='OLD' ) 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/cfg/variable.cfg

    r6140 r7339  
    11# name       | units          | axis | pt| interpolation   | long name                             | standard name                                   
    2 X            | 1              | X    |   |                 |                                       | projection_x_coordinate                
    3 Y            | 1              | Y    |   |                 |                                       | projection_y_coordinate                
    4 Z            | 1              | Z    |   |                 |                                       | projection_z_coordinate                
    5 T            | 1              | T    |   |                 |                                       | projection_t_coordinate                
     2X            | unitless       | X    |   |                 |                                       | projection_x_coordinate                
     3Y            | unitless       | Y    |   |                 |                                       | projection_y_coordinate                
     4Z            | unitless       | Z    |   |                 |                                       | projection_z_coordinate                
     5T            | unitless       | T    |   |                 |                                       | projection_t_coordinate                
    66nav_lon      | degrees_east   | XY   | T | cubic           | Longitude                             | longitude                                    
    77nav_lat      | degrees_north  | XY   | T | cubic           | Latitude                              | latitude                          
     
    4343kt           |                |      |   |                 |                                       |                                   
    4444rdt          |                |      |   |                 |                                       |                                   
     45rdttra1      |                |      |   |                 |                                       |                                   
    4546utau_b       |                | XY   | U |                 |                                       |surface_downward_eastward_stress   
    4647vtau_b       |                | XY   | V |                 |                                       |surface_downward_northward_stress  
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/Doxyfile

    r5037 r7339  
    4545# quick idea about the purpose of the project. Keep the description short. 
    4646 
    47 PROJECT_BRIEF          = "System and Interface for oceanic RElocable Nesting" 
     47PROJECT_BRIEF          = "System and Interface for oceanic RElocatable Nesting" 
    4848 
    4949# With the PROJECT_LOGO tag one can specify an logo or icon that is included in 
     
    20692069# The default value is: NO. 
    20702070 
    2071 HAVE_DOT               = YES 
     2071HAVE_DOT               = NO 
    20722072 
    20732073# The DOT_NUM_THREADS specifies the number of dot invocations doxygen is allowed 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/attribute.f90

    r5617 r7339  
    8383!> @date November, 2014  
    8484!> - Fix memory leaks bug 
     85!> @date September, 2015 
     86!> - manage useless (dummy) attributes 
    8587! 
    8688!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    98100   PUBLIC :: TATT       !< attribute structure 
    99101 
     102   PRIVATE :: cm_dumatt !< dummy attribute array 
     103 
    100104   ! function and subroutine 
    101105   PUBLIC :: att_init       !< initialize attribute structure 
     
    105109   PUBLIC :: att_get_index  !< get attribute index, in an array of attribute structure 
    106110   PUBLIC :: att_get_id     !< get attribute id, read from file 
     111   PUBLIC :: att_get_dummy  !< fill dummy attribute array 
     112   PUBLIC :: att_is_dummy   !< check if attribute is defined as dummy attribute 
    107113 
    108114   PRIVATE :: att__clean_unit ! clean attribute strcuture 
     
    135141   END TYPE TATT 
    136142 
     143   CHARACTER(LEN=lc), DIMENSION(ip_maxdum), SAVE :: cm_dumatt !< dummy attribute 
     144 
    137145   INTERFACE att_init 
    138146      MODULE PROCEDURE att__init_c     
     
    12511259 
    12521260   END SUBROUTINE att__clean_arr 
     1261   !------------------------------------------------------------------- 
     1262   !> @brief This subroutine fill dummy attribute array 
     1263   ! 
     1264   !> @author J.Paul 
     1265   !> @date September, 2015 - Initial Version 
     1266   !> @date Marsh, 2016 
     1267   !> - close file (bugfix) 
     1268   ! 
     1269   !> @param[in] cd_dummy dummy configuration file 
     1270   !------------------------------------------------------------------- 
     1271   SUBROUTINE att_get_dummy( cd_dummy ) 
     1272      IMPLICIT NONE 
     1273      ! Argument 
     1274      CHARACTER(LEN=*), INTENT(IN) :: cd_dummy 
     1275 
     1276      ! local variable 
     1277      INTEGER(i4)   :: il_fileid 
     1278      INTEGER(i4)   :: il_status 
     1279 
     1280      LOGICAL       :: ll_exist 
     1281 
     1282      ! loop indices 
     1283      ! namelist 
     1284      CHARACTER(LEN=lc), DIMENSION(ip_maxdum) :: cn_dumvar 
     1285      CHARACTER(LEN=lc), DIMENSION(ip_maxdum) :: cn_dumdim 
     1286      CHARACTER(LEN=lc), DIMENSION(ip_maxdum) :: cn_dumatt 
     1287 
     1288      !---------------------------------------------------------------- 
     1289      NAMELIST /namdum/ &   !< dummy namelist 
     1290      &  cn_dumvar, &       !< variable  name 
     1291      &  cn_dumdim, &       !< dimension name 
     1292      &  cn_dumatt          !< attribute name 
     1293      !---------------------------------------------------------------- 
     1294 
     1295      ! init 
     1296      cm_dumatt(:)='' 
     1297 
     1298      ! read namelist 
     1299      INQUIRE(FILE=TRIM(cd_dummy), EXIST=ll_exist) 
     1300      IF( ll_exist )THEN 
     1301     
     1302         il_fileid=fct_getunit() 
     1303    
     1304         OPEN( il_fileid, FILE=TRIM(cd_dummy), & 
     1305         &                FORM='FORMATTED',       & 
     1306         &                ACCESS='SEQUENTIAL',    & 
     1307         &                STATUS='OLD',           & 
     1308         &                ACTION='READ',          & 
     1309         &                IOSTAT=il_status) 
     1310         CALL fct_err(il_status) 
     1311         IF( il_status /= 0 )THEN 
     1312            CALL logger_fatal("DIM GET DUMMY: opening "//TRIM(cd_dummy)) 
     1313         ENDIF 
     1314    
     1315         READ( il_fileid, NML = namdum ) 
     1316         cm_dumatt(:)=cn_dumatt(:) 
     1317 
     1318         CLOSE( il_fileid ) 
     1319 
     1320      ENDIF 
     1321    
     1322   END SUBROUTINE att_get_dummy 
     1323   !------------------------------------------------------------------- 
     1324   !> @brief This function check if attribute is defined as dummy attribute 
     1325   !> in configuraton file 
     1326   !> 
     1327   !> @author J.Paul 
     1328   !> @date September, 2015 - Initial Version 
     1329   ! 
     1330   !> @param[in] td_att attribute structure 
     1331   !> @return true if attribute is dummy attribute 
     1332   !------------------------------------------------------------------- 
     1333   FUNCTION att_is_dummy(td_att) 
     1334      IMPLICIT NONE 
     1335 
     1336      ! Argument       
     1337      TYPE(TATT), INTENT(IN) :: td_att 
     1338       
     1339      ! function 
     1340      LOGICAL :: att_is_dummy 
     1341       
     1342      ! loop indices 
     1343      INTEGER(i4) :: ji 
     1344      !---------------------------------------------------------------- 
     1345 
     1346      att_is_dummy=.FALSE. 
     1347      DO ji=1,ip_maxdum 
     1348         IF( fct_lower(td_att%c_name) == fct_lower(cm_dumatt(ji)) )THEN 
     1349            att_is_dummy=.TRUE. 
     1350            EXIT 
     1351         ENDIF 
     1352      ENDDO 
     1353 
     1354   END FUNCTION att_is_dummy 
    12531355END MODULE att 
    12541356 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/boundary.f90

    r5609 r7339  
    482482   !> ex : cn_north='index1,first1,last1(width)|index2,first2,last2' 
    483483   !> 
    484    !> @note Boundaries are compute on T point, but expressed on U,V point. 
     484   !> @warn Boundaries are compute on T point, but expressed on U,V point. 
    485485   !> change will be done to get data on other point when need be.  
    486486   !> 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/create_bathy.f90

    r5617 r7339  
    88!> @file 
    99!> @brief  
    10 !> This program create fine grid bathymetry file. 
     10!> This program creates fine grid bathymetry file. 
    1111!> 
    1212!> @details 
     
    2727!>    you could find a template of the namelist in templates directory. 
    2828!> 
    29 !>    create_bathy.nam comprise 7 namelists:<br/> 
     29!>    create_bathy.nam contains 7 namelists:<br/> 
    3030!>       - logger namelist (namlog) 
    3131!>       - config namelist (namcfg) 
     
    3636!>       - output namelist (namout) 
    3737!>     
    38 !>    @note  
    39 !>       All namelists have to be in file create_bathy.nam, however variables of 
    40 !>       those namelists are all optional. 
    41 !> 
    4238!>    * _logger namelist (namlog)_:<br/> 
    4339!>       - cn_logfile   : log filename 
     
    4945!>       - cn_varcfg : variable configuration file  
    5046!> (see ./SIREN/cfg/variable.cfg) 
     47!>       - cn_dumcfg : useless (dummy) configuration file, for useless  
     48!> dimension or variable (see ./SIREN/cfg/dummy.cfg). 
    5149!> 
    5250!>    * _coarse grid namelist (namcrs)_:<br/> 
     
    6159!> 
    6260!>    * _variable namelist (namvar)_:<br/> 
    63 !>       - cn_varinfo : list of variable and extra information about request(s)  
    64 !>       to be used.<br/> 
    65 !>          each elements of *cn_varinfo* is a string character 
    66 !>          (separated by ',').<br/> 
    67 !>          it is composed of the variable name follow by ':',  
    68 !>          then request(s) to be used on this variable.<br/>  
    69 !>          request could be: 
    70 !>             - int = interpolation method 
    71 !>             - ext = extrapolation method 
    72 !>             - flt = filter method 
    73 !>             - min = minimum value 
    74 !>             - max = maximum value 
    75 !>             - unt = new units 
    76 !>             - unf = unit scale factor (linked to new units) 
    77 !> 
    78 !>                requests must be separated by ';'.<br/> 
    79 !>                order of requests does not matter.<br/> 
    80 !> 
    81 !>          informations about available method could be find in @ref interp, 
    82 !>          @ref extrap and @ref filter modules.<br/> 
    83 !>          Example: 'Bathymetry: flt=2*hamming(2,3); min=0' 
    84 !>          @note  
    85 !>             If you do not specify a method which is required,  
    86 !>             default one is apply. 
    87 !>          @warning  
    88 !>             variable name must be __Bathymetry__ here. 
    8961!>       - cn_varfile : list of variable, and corresponding file.<br/>  
    9062!>          *cn_varfile* is the path and filename of the file where find 
     
    10880!>             - 'Bathymetry:5000,5000,5000/5000,3000,5000/5000,5000,5000' 
    10981!> 
     82!>       - cn_varinfo : list of variable and extra information about request(s)  
     83!>       to be used.<br/> 
     84!>          each elements of *cn_varinfo* is a string character 
     85!>          (separated by ',').<br/> 
     86!>          it is composed of the variable name follow by ':',  
     87!>          then request(s) to be used on this variable.<br/>  
     88!>          request could be: 
     89!>             - int = interpolation method 
     90!>             - ext = extrapolation method 
     91!>             - flt = filter method 
     92!>             - min = minimum value 
     93!>             - max = maximum value 
     94!>             - unt = new units 
     95!>             - unf = unit scale factor (linked to new units) 
     96!> 
     97!>                requests must be separated by ';'.<br/> 
     98!>                order of requests does not matter.<br/> 
     99!> 
     100!>          informations about available method could be find in @ref interp, 
     101!>          @ref extrap and @ref filter modules.<br/> 
     102!>          Example: 'Bathymetry: flt=2*hamming(2,3); min=0' 
     103!>          @note  
     104!>             If you do not specify a method which is required,  
     105!>             default one is apply. 
     106!>          @warning  
     107!>             variable name must be __Bathymetry__ here. 
     108!> 
    110109!>    * _nesting namelist (namnst)_:<br/> 
    111110!>       - in_rhoi  : refinement factor in i-direction 
     
    127126!> - extrapolate all land points. 
    128127!> - allow to change unit. 
     128!> @date September, 2015 
     129!> - manage useless (dummy) variable, attributes, and dimension 
     130!> @date January,2016 
     131!> - add create_bathy_check_depth as in create_boundary 
     132!> - add create_bathy_check_time  as in create_boundary 
     133!> @date February, 2016 
     134!> - do not closed sea for east-west cyclic domain 
    129135! 
    130136!> @todo 
    131 !> - use create_bathy_check_depth as in create_boundary 
    132 !> - use create_bathy_check_time  as in create_boundary 
    133137!> - check tl_multi is not empty 
    134138!> 
     
    167171   INTEGER(i4)                                        :: il_status 
    168172   INTEGER(i4)                                        :: il_fileid 
    169    INTEGER(i4)                                        :: il_varid 
    170173   INTEGER(i4)                                        :: il_attid 
    171174   INTEGER(i4)                                        :: il_imin0 
     
    179182 
    180183   LOGICAL                                            :: ll_exist 
     184   LOGICAL                                            :: ll_fillclosed 
    181185 
    182186   TYPE(TMPP)                                         :: tl_coord0 
     
    208212   ! namelist variable 
    209213   ! namlog 
    210    CHARACTER(LEN=lc) :: cn_logfile = 'create_bathy.log'  
    211    CHARACTER(LEN=lc) :: cn_verbosity = 'warning'  
    212    INTEGER(i4)       :: in_maxerror = 5 
     214   CHARACTER(LEN=lc)                       :: cn_logfile = 'create_bathy.log'  
     215   CHARACTER(LEN=lc)                       :: cn_verbosity = 'warning'  
     216   INTEGER(i4)                             :: in_maxerror = 5 
    213217 
    214218   ! namcfg 
    215    CHARACTER(LEN=lc) :: cn_varcfg = 'variable.cfg'  
     219   CHARACTER(LEN=lc)                       :: cn_varcfg = './cfg/variable.cfg'  
     220   CHARACTER(LEN=lc)                       :: cn_dumcfg = './cfg/dummy.cfg'  
    216221 
    217222   ! namcrs 
    218    CHARACTER(LEN=lc) :: cn_coord0 = ''  
    219    INTEGER(i4)       :: in_perio0 = -1 
     223   CHARACTER(LEN=lc)                       :: cn_coord0 = ''  
     224   INTEGER(i4)                             :: in_perio0 = -1 
    220225 
    221226   ! namfin 
    222    CHARACTER(LEN=lc) :: cn_coord1 = '' 
    223    INTEGER(i4)       :: in_perio1 = -1 
    224    LOGICAL           :: ln_fillclosed = .TRUE. 
     227   CHARACTER(LEN=lc)                       :: cn_coord1 = '' 
     228   INTEGER(i4)                             :: in_perio1 = -1 
     229   LOGICAL                                 :: ln_fillclosed = .TRUE. 
    225230 
    226231   ! namvar 
     232   CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varfile = '' 
    227233   CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varinfo = '' 
    228    CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varfile = '' 
    229234 
    230235   ! namnst 
    231    INTEGER(i4)       :: in_rhoi  = 1 
    232    INTEGER(i4)       :: in_rhoj  = 1 
     236   INTEGER(i4)                             :: in_rhoi  = 1 
     237   INTEGER(i4)                             :: in_rhoj  = 1 
    233238 
    234239   ! namout 
    235    CHARACTER(LEN=lc) :: cn_fileout = 'bathy_fine.nc'  
     240   CHARACTER(LEN=lc)                       :: cn_fileout = 'bathy_fine.nc'  
    236241   !------------------------------------------------------------------- 
    237242 
     
    242247 
    243248   NAMELIST /namcfg/ &   !< configuration namelist 
    244    &  cn_varcfg          !< variable configuration file 
     249   &  cn_varcfg, &       !< variable configuration file 
     250   &  cn_dumcfg          !< dummy configuration file 
    245251 
    246252   NAMELIST /namcrs/ &   !< coarse grid namelist 
     
    254260  
    255261   NAMELIST /namvar/ &   !< variable namelist 
    256    &  cn_varinfo, &      !< list of variable and interpolation method to be used. (ex: 'votemper:linear','vosaline:cubic' ) 
    257    &  cn_varfile         !< list of variable file 
     262   &  cn_varfile, &      !< list of variable file 
     263   &  cn_varinfo         !< list of variable and interpolation method to be used. (ex: 'votemper:linear','vosaline:cubic' ) 
    258264    
    259265   NAMELIST /namnst/ &   !< nesting namelist 
     
    302308      CALL var_def_extra(TRIM(cn_varcfg)) 
    303309 
     310      ! get dummy variable 
     311      CALL var_get_dummy(TRIM(cn_dumcfg)) 
     312      ! get dummy dimension 
     313      CALL dim_get_dummy(TRIM(cn_dumcfg)) 
     314      ! get dummy attribute 
     315      CALL att_get_dummy(TRIM(cn_dumcfg)) 
     316 
    304317      READ( il_fileid, NML = namcrs ) 
    305318      READ( il_fileid, NML = namfin ) 
     
    309322      ! match variable with file 
    310323      tl_multi=multi_init(cn_varfile) 
    311        
     324  
    312325      READ( il_fileid, NML = namnst ) 
    313326      READ( il_fileid, NML = namout ) 
     
    322335 
    323336      PRINT *,"ERROR in create_bathy: can't find "//TRIM(cl_namelist) 
     337      STOP 
    324338 
    325339   ENDIF 
     
    343357      &     "check namelist") 
    344358   ENDIF 
     359 
     360   ! do not closed sea for east-west cyclic domain 
     361   ll_fillclosed=ln_fillclosed 
     362   IF( tl_coord1%i_perio == 1 ) ll_fillclosed=.FALSE. 
    345363 
    346364   ! check 
     
    417435 
    418436            ! get or check depth value 
    419             IF( tl_multi%t_mpp(ji)%t_proc(1)%i_depthid /= 0 )THEN 
    420                il_varid=tl_multi%t_mpp(ji)%t_proc(1)%i_depthid 
    421                IF( ASSOCIATED(tl_depth%d_value) )THEN 
    422                   tl_tmp=iom_mpp_read_var(tl_mpp,il_varid) 
    423                   IF( ANY( tl_depth%d_value(:,:,:,:) /= & 
    424                   &        tl_tmp%d_value(:,:,:,:) ) )THEN 
    425                      CALL logger_fatal("CREATE BATHY: depth value from "//& 
    426                      &  TRIM(tl_multi%t_mpp(ji)%c_name)//" not conform "//& 
    427                      &  " to those from former file(s).") 
    428                   ENDIF 
    429                   CALL var_clean(tl_tmp) 
    430                ELSE 
    431                   tl_depth=iom_mpp_read_var(tl_mpp,il_varid) 
    432                ENDIF 
    433             ENDIF 
     437            CALL create_bathy_check_depth( tl_mpp, tl_depth ) 
    434438 
    435439            ! get or check time value 
    436             IF( tl_multi%t_mpp(ji)%t_proc(1)%i_timeid /= 0 )THEN 
    437                il_varid=tl_multi%t_mpp(ji)%t_proc(1)%i_timeid 
    438                IF( ASSOCIATED(tl_time%d_value) )THEN 
    439                   tl_tmp=iom_mpp_read_var(tl_mpp,il_varid) 
    440                   IF( ANY( tl_time%d_value(:,:,:,:) /= & 
    441                   &        tl_tmp%d_value(:,:,:,:) ) )THEN 
    442                      CALL logger_fatal("CREATE BATHY: time value from "//& 
    443                      &  TRIM(tl_multi%t_mpp(ji)%c_name)//" not conform "//& 
    444                      &  " to those from former file(s).") 
    445                   ENDIF 
    446                   CALL var_clean(tl_tmp) 
    447                ELSE 
    448                   tl_time=iom_mpp_read_var(tl_mpp,il_varid) 
    449                ENDIF 
    450             ENDIF 
     440            CALL create_bathy_check_time( tl_mpp, tl_time ) 
    451441 
    452442            ! close mpp file 
    453443            CALL iom_mpp_close(tl_mpp) 
    454444 
    455             IF( ANY( tl_mpp%t_dim(1:2)%i_len /= & 
    456             &        tl_coord0%t_dim(1:2)%i_len) )THEN 
     445            IF( ANY(tl_mpp%t_dim(1:2)%i_len /= tl_coord0%t_dim(1:2)%i_len).OR.& 
     446            &   ALL(il_rho(:)==1) )THEN 
    457447               !- extract bathymetry from fine grid bathymetry  
    458448               DO jj=1,tl_multi%t_mpp(ji)%t_proc(1)%i_nvar 
     
    505495 
    506496         ! fill closed sea 
    507          IF( ln_fillclosed )THEN 
     497         IF( ll_fillclosed )THEN 
    508498            ALLOCATE( il_mask(tl_var(jk)%t_dim(1)%i_len, & 
    509499            &                 tl_var(jk)%t_dim(2)%i_len) ) 
     
    526516         &   dl_minbat <= 0._dp  )THEN 
    527517            CALL logger_debug("CREATE BATHY: min value "//TRIM(fct_str(dl_minbat))) 
    528             CALL logger_error("CREATE BATHY: Bathymetry has value <= 0") 
     518            CALL logger_fatal("CREATE BATHY: Bathymetry has value <= 0") 
    529519         ENDIF 
    530520 
     
    973963      CALL dom_del_extra( tl_var, tl_dom, il_rho(:) ) 
    974964 
     965      CALL dom_clean_extra( tl_dom ) 
     966 
    975967      !- add ghost cell 
    976968      CALL grid_add_ghost(tl_var,tl_dom%i_ghost(:,:)) 
     
    11091101 
    11101102   END SUBROUTINE create_bathy_interp 
     1103   !------------------------------------------------------------------- 
     1104   !> @brief 
     1105   !> This subroutine get depth variable value in an open mpp structure 
     1106   !> and check if agree with already input depth variable. 
     1107   !>  
     1108   !> @details  
     1109   !> 
     1110   !> @author J.Paul 
     1111   !> @date January, 2016 - Initial Version 
     1112   !> 
     1113   !> @param[in] td_mpp       mpp structure 
     1114   !> @param[inout] td_depth  depth variable structure  
     1115   !------------------------------------------------------------------- 
     1116   SUBROUTINE create_bathy_check_depth( td_mpp, td_depth ) 
     1117 
     1118      IMPLICIT NONE 
     1119 
     1120      ! Argument 
     1121      TYPE(TMPP) , INTENT(IN   ) :: td_mpp 
     1122      TYPE(TVAR) , INTENT(INOUT) :: td_depth 
     1123 
     1124      ! local variable 
     1125      INTEGER(i4) :: il_varid 
     1126      TYPE(TVAR)  :: tl_depth 
     1127      ! loop indices 
     1128      !---------------------------------------------------------------- 
     1129 
     1130      ! get or check depth value 
     1131      IF( td_mpp%t_proc(1)%i_depthid /= 0 )THEN 
     1132 
     1133         il_varid=td_mpp%t_proc(1)%i_depthid 
     1134         IF( ASSOCIATED(td_depth%d_value) )THEN 
     1135 
     1136            tl_depth=iom_mpp_read_var(td_mpp, il_varid) 
     1137 
     1138            IF( ANY( td_depth%d_value(:,:,:,:) /= & 
     1139            &        tl_depth%d_value(:,:,:,:) ) )THEN 
     1140 
     1141               CALL logger_warn("CREATE BATHY: depth value from "//& 
     1142               &  TRIM(td_mpp%c_name)//" not conform "//& 
     1143               &  " to those from former file(s).") 
     1144 
     1145            ENDIF 
     1146            CALL var_clean(tl_depth) 
     1147 
     1148         ELSE 
     1149            td_depth=iom_mpp_read_var(td_mpp,il_varid) 
     1150         ENDIF 
     1151 
     1152      ENDIF 
     1153       
     1154   END SUBROUTINE create_bathy_check_depth 
     1155   !------------------------------------------------------------------- 
     1156   !> @brief 
     1157   !> This subroutine get date and time in an open mpp structure 
     1158   !> and check if agree with date and time already read. 
     1159   !>  
     1160   !> @details  
     1161   !> 
     1162   !> @author J.Paul 
     1163   !> @date January, 2016 - Initial Version 
     1164   !> 
     1165   !> @param[in] td_mpp      mpp structure 
     1166   !> @param[inout] td_time  time variable structure  
     1167   !------------------------------------------------------------------- 
     1168   SUBROUTINE create_bathy_check_time( td_mpp, td_time ) 
     1169 
     1170      IMPLICIT NONE 
     1171 
     1172      ! Argument 
     1173      TYPE(TMPP), INTENT(IN   ) :: td_mpp 
     1174      TYPE(TVAR), INTENT(INOUT) :: td_time 
     1175 
     1176      ! local variable 
     1177      INTEGER(i4) :: il_varid 
     1178      TYPE(TVAR)  :: tl_time 
     1179 
     1180      TYPE(TDATE) :: tl_date1 
     1181      TYPE(TDATE) :: tl_date2 
     1182      ! loop indices 
     1183      !---------------------------------------------------------------- 
     1184 
     1185      ! get or check depth value 
     1186      IF( td_mpp%t_proc(1)%i_timeid /= 0 )THEN 
     1187 
     1188         il_varid=td_mpp%t_proc(1)%i_timeid 
     1189         IF( ASSOCIATED(td_time%d_value) )THEN 
     1190 
     1191            tl_time=iom_mpp_read_var(td_mpp, il_varid) 
     1192 
     1193            tl_date1=var_to_date(td_time) 
     1194            tl_date2=var_to_date(tl_time) 
     1195            IF( tl_date1 - tl_date2 /= 0 )THEN 
     1196 
     1197               CALL logger_warn("CREATE BATHY: date from "//& 
     1198               &  TRIM(td_mpp%c_name)//" not conform "//& 
     1199               &  " to those from former file(s).") 
     1200 
     1201            ENDIF 
     1202            CALL var_clean(tl_time) 
     1203 
     1204         ELSE 
     1205            td_time=iom_mpp_read_var(td_mpp,il_varid) 
     1206         ENDIF 
     1207 
     1208      ENDIF 
     1209       
     1210   END SUBROUTINE create_bathy_check_time 
    11111211END PROGRAM create_bathy 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/create_coord.f90

    r5609 r7339  
    99!> @file 
    1010!> @brief  
    11 !> This program create fine grid coordinate file. 
     11!> This program creates fine grid coordinate file. 
    1212!> 
    1313!> @details 
     
    2727!>    you could find a template of the namelist in templates directory. 
    2828!> 
    29 !>    create_coord.nam comprise 6 namelists:<br/> 
     29!>    create_coord.nam contains 6 namelists:<br/> 
    3030!>       - logger namelist (namlog) 
    3131!>       - config namelist (namcfg) 
     
    3535!>       - output namelist (namout) 
    3636!>     
    37 !>    @note  
    38 !>       All namelists have to be in file create_coord.nam,  
    39 !>       however variables of those namelists are all optional. 
    40 !> 
    4137!>    * _logger namelist (namlog)_:<br/> 
    4238!>       - cn_logfile   : log filename 
     
    4844!>       - cn_varcfg : variable configuration file  
    4945!> (see ./SIREN/cfg/variable.cfg) 
     46!>       - cn_dumcfg : useless (dummy) configuration file, for useless  
     47!> dimension or variable (see ./SIREN/cfg/dummy.cfg). 
    5048!> 
    5149!>    * _coarse grid namelist (namcrs)_:<br/> 
     
    6462!>             - int = interpolation method 
    6563!>             - ext = extrapolation method 
    66 !>             - flt = filter method 
    6764!>  
    6865!>                requests must be separated by ';' .<br/> 
     
    7269!>          @ref extrap and @ref filter modules.<br/> 
    7370!> 
    74 !>          Example: 'votemper: int=linear; flt=hann(2,3); ext=dist_weight',  
    75 !>          'vosaline: int=cubic'<br/> 
     71!>          Example: 'glamt: int=linear; ext=dist_weight',  
     72!>          'e1t: int=cubic/rhoi'<br/> 
    7673!>          @note  
    7774!>             If you do not specify a method which is required,  
     
    103100!> - compute offset considering grid point 
    104101!> - add global attributes in output file 
     102!> @date September, 2015 
     103!> - manage useless (dummy) variable, attributes, and dimension 
    105104!> 
    106105!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    167166 
    168167   ! namcfg 
    169    CHARACTER(LEN=lc) :: cn_varcfg = '../cfg/variable.cfg'  
     168   CHARACTER(LEN=lc) :: cn_varcfg = './cfg/variable.cfg'  
     169   CHARACTER(LEN=lc) :: cn_dumcfg = './cfg/dummy.cfg' 
    170170 
    171171   ! namcrs 
     
    194194 
    195195   NAMELIST /namcfg/ &  !  config namelist 
    196    &  cn_varcfg         !< variable configuration file 
     196   &  cn_varcfg, &       !< variable configuration file 
     197   &  cn_dumcfg          !< dummy configuration file 
    197198 
    198199   NAMELIST /namcrs/ &  !  coarse grid namelist 
     
    254255      CALL var_def_extra(TRIM(cn_varcfg)) 
    255256 
     257      ! get dummy variable 
     258      CALL var_get_dummy(TRIM(cn_dumcfg)) 
     259      ! get dummy dimension 
     260      CALL dim_get_dummy(TRIM(cn_dumcfg)) 
     261      ! get dummy attribute 
     262      CALL att_get_dummy(TRIM(cn_dumcfg)) 
     263 
    256264      READ( il_fileid, NML = namcrs ) 
    257265      READ( il_fileid, NML = namvar ) 
     
    354362   ENDDO 
    355363 
     364   ! clean 
     365   CALL dom_clean_extra( tl_dom ) 
     366 
    356367   ! close mpp files 
    357368   CALL iom_dom_close(tl_coord0) 
     
    388399   CALL file_add_att(tl_fileout, tl_att)    
    389400 
    390    tl_att=att_init("src_i_indices",(/in_imin0,in_imax0/)) 
     401   tl_att=att_init("src_i_indices",(/tl_dom%i_imin,tl_dom%i_imax/)) 
    391402   CALL file_add_att(tl_fileout, tl_att)    
    392    tl_att=att_init("src_j_indices",(/in_jmin0,in_jmax0/)) 
     403   tl_att=att_init("src_j_indices",(/tl_dom%i_jmin,tl_dom%i_jmax/)) 
    393404   CALL file_add_att(tl_fileout, tl_att) 
    394405   IF( .NOT. ALL(il_rho(:)==1) )THEN 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/create_restart.f90

    r5617 r7339  
    99!> @file 
    1010!> @brief  
    11 !> This program create restart file. 
     11!> This program creates restart file. 
    1212!> 
    1313!> @details 
    1414!> @section sec1 method 
    1515!> Variables could be extracted from fine grid file, interpolated from coarse 
    16 !> grid file or restart file, or manually written.<br/>  
    17 !> Then they are split over new decomposition.  
     16!> grid file or restart file. Variables could also be manually written.<br/>  
     17!> Then they are split over new layout.  
    1818!> @note  
    1919!>    method could be different for each variable. 
     
    2828!>    you could find a template of the namelist in templates directory. 
    2929!> 
    30 !>    create_restart.nam comprise 9 namelists:<br/> 
     30!>    create_restart.nam contains 9 namelists:<br/> 
    3131!>       - logger namelist (namlog) 
    3232!>       - config namelist (namcfg) 
     
    3939!>       - output namelist (namout) 
    4040!>     
    41 !>    @note  
    42 !>       All namelists have to be in file create_restart.nam  
    43 !>       however variables of those namelists are all optional. 
    44 !> 
    4541!>    * _logger namelist (namlog)_:<br/> 
    4642!>       - cn_logfile   : log filename 
     
    5248!>       - cn_varcfg : variable configuration file 
    5349!> (see ./SIREN/cfg/variable.cfg) 
     50!>       - cn_dumcfg : useless (dummy) configuration file, for useless  
     51!> dimension or variable (see ./SIREN/cfg/dummy.cfg). 
    5452!> 
    5553!>    * _coarse grid namelist (namcrs):<br/> 
     
    8280!> 
    8381!>    * _variable namelist (namvar)_:<br/> 
    84 !>       - cn_varinfo : list of variable and extra information about request(s)  
    85 !>       to be used.<br/> 
    86 !>          each elements of *cn_varinfo* is a string character 
    87 !>          (separated by ',').<br/> 
    88 !>          it is composed of the variable name follow by ':',  
    89 !>          then request(s) to be used on this variable.<br/>  
    90 !>          request could be: 
    91 !>             - int = interpolation method 
    92 !>             - ext = extrapolation method 
    93 !>             - flt = filter method 
    94 !>             - min = minimum value 
    95 !>             - max = maximum value 
    96 !>             - unt = new units 
    97 !>             - unf = unit scale factor (linked to new units) 
    98 !> 
    99 !>             requests must be separated by ';'.<br/> 
    100 !>             order of requests does not matter.<br/> 
    101 !> 
    102 !>          informations about available method could be find in @ref interp, 
    103 !>          @ref extrap and @ref filter.<br/> 
    104 !>          Example: 'votemper: int=linear; flt=hann; ext=dist_weight','vosaline: int=cubic' 
    105 !>          @note  
    106 !>             If you do not specify a method which is required,  
    107 !>             default one is apply. 
    108 !>       - cn_varfile : list of variable, and corresponding file<br/>  
     82!>       - cn_varfile : list of variable, and associated file<br/>  
    10983!>          *cn_varfile* is the path and filename of the file where find 
    11084!>          variable.<br/> 
     
    131105!>             - 'all:restart.dimg' 
    132106!> 
     107!>       - cn_varinfo : list of variable and extra information about request(s)  
     108!>       to be used.<br/> 
     109!>          each elements of *cn_varinfo* is a string character 
     110!>          (separated by ',').<br/> 
     111!>          it is composed of the variable name follow by ':',  
     112!>          then request(s) to be used on this variable.<br/>  
     113!>          request could be: 
     114!>             - int = interpolation method 
     115!>             - ext = extrapolation method 
     116!>             - flt = filter method 
     117!>             - min = minimum value 
     118!>             - max = maximum value 
     119!>             - unt = new units 
     120!>             - unf = unit scale factor (linked to new units) 
     121!> 
     122!>             requests must be separated by ';'.<br/> 
     123!>             order of requests does not matter.<br/> 
     124!> 
     125!>          informations about available method could be find in @ref interp, 
     126!>          @ref extrap and @ref filter.<br/> 
     127!>          Example: 'votemper: int=linear; flt=hann; ext=dist_weight', 
     128!>                   'vosaline: int=cubic' 
     129!>          @note  
     130!>             If you do not specify a method which is required,  
     131!>             default one is apply. 
     132!> 
    133133!>    * _nesting namelist (namnst)_:<br/> 
    134134!>       - in_rhoi  : refinement factor in i-direction 
    135135!>       - in_rhoj  : refinement factor in j-direction 
    136136!>       @note  
    137 !>          coarse grid indices will be deduced from fine grid 
     137!>          coarse grid indices will be computed from fine grid 
    138138!>          coordinate file. 
    139139!> 
     
    141141!>       - cn_fileout : output file 
    142142!>       - ln_extrap : extrapolate land point or not 
    143 !>       - in_niproc : i-direction number of processor 
    144 !>       - in_njproc : j-direction numebr of processor 
     143!>       - in_niproc : number of processor in i-direction 
     144!>       - in_njproc : number of processor in j-direction 
    145145!>       - in_nproc  : total number of processor to be used 
    146146!>       - cn_type   : output format ('dimg', 'cdf') 
     
    156156!> - extrapolate all land points, and add ln_extrap in namelist. 
    157157!> - allow to change unit. 
     158!> @date September, 2015 
     159!> - manage useless (dummy) variable, attributes, and dimension 
    158160!> 
    159161!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    173175   USE iom                             ! I/O manager 
    174176   USE grid                            ! grid manager 
    175    USE vgrid                            ! vertical grid manager 
     177   USE vgrid                           ! vertical grid manager 
    176178   USE extrap                          ! extrapolation manager 
    177179   USE interp                          ! interpolation manager 
     
    183185 
    184186   IMPLICIT NONE 
    185  
    186187 
    187188   ! local variable 
     
    212213 
    213214   LOGICAL                                            :: ll_exist 
     215   LOGICAL                                            :: ll_sameGrid 
    214216 
    215217   TYPE(TDOM)                                         :: tl_dom1 
     
    242244   ! namelist variable 
    243245   ! namlog 
    244    CHARACTER(LEN=lc) :: cn_logfile = 'create_restart.log'  
    245    CHARACTER(LEN=lc) :: cn_verbosity = 'warning'  
    246    INTEGER(i4)       :: in_maxerror = 5 
     246   CHARACTER(LEN=lc)                       :: cn_logfile = 'create_restart.log'  
     247   CHARACTER(LEN=lc)                       :: cn_verbosity = 'warning'  
     248   INTEGER(i4)                             :: in_maxerror = 5 
    247249 
    248250   ! namcfg 
    249    CHARACTER(LEN=lc) :: cn_varcfg = 'variable.cfg'  
     251   CHARACTER(LEN=lc)                       :: cn_varcfg = 'variable.cfg'  
     252   CHARACTER(LEN=lc)                       :: cn_dumcfg = 'dummy.cfg' 
    250253 
    251254   ! namcrs 
    252    CHARACTER(LEN=lc) :: cn_coord0 = ''  
    253    INTEGER(i4)       :: in_perio0 = -1 
     255   CHARACTER(LEN=lc)                       :: cn_coord0 = ''  
     256   INTEGER(i4)                             :: in_perio0 = -1 
    254257 
    255258   ! namfin 
    256    CHARACTER(LEN=lc) :: cn_coord1 = '' 
    257    CHARACTER(LEN=lc) :: cn_bathy1 = '' 
    258    INTEGER(i4)       :: in_perio1 = -1 
     259   CHARACTER(LEN=lc)                       :: cn_coord1 = '' 
     260   CHARACTER(LEN=lc)                       :: cn_bathy1 = '' 
     261   INTEGER(i4)                             :: in_perio1 = -1 
    259262 
    260263   !namzgr 
    261    REAL(dp)          :: dn_pp_to_be_computed = 0._dp 
    262    REAL(dp)          :: dn_ppsur     = -3958.951371276829_dp 
    263    REAL(dp)          :: dn_ppa0      =   103.9530096000000_dp 
    264    REAL(dp)          :: dn_ppa1      =     2.4159512690000_dp 
    265    REAL(dp)          :: dn_ppa2      =   100.7609285000000_dp 
    266    REAL(dp)          :: dn_ppkth     =    15.3510137000000_dp 
    267    REAL(dp)          :: dn_ppkth2    =    48.0298937200000_dp 
    268    REAL(dp)          :: dn_ppacr     =     7.0000000000000_dp 
    269    REAL(dp)          :: dn_ppacr2    =    13.000000000000_dp 
    270    REAL(dp)          :: dn_ppdzmin  = 6._dp 
    271    REAL(dp)          :: dn_pphmax    = 5750._dp 
    272    INTEGER(i4)       :: in_nlevel    = 75 
     264   REAL(dp)                                :: dn_pp_to_be_computed = 0._dp 
     265   REAL(dp)                                :: dn_ppsur   = -3958.951371276829_dp 
     266   REAL(dp)                                :: dn_ppa0    =   103.953009600000_dp 
     267   REAL(dp)                                :: dn_ppa1    =     2.415951269000_dp 
     268   REAL(dp)                                :: dn_ppa2    =   100.760928500000_dp 
     269   REAL(dp)                                :: dn_ppkth   =    15.351013700000_dp 
     270   REAL(dp)                                :: dn_ppkth2  =    48.029893720000_dp 
     271   REAL(dp)                                :: dn_ppacr   =     7.000000000000_dp 
     272   REAL(dp)                                :: dn_ppacr2  =    13.000000000000_dp 
     273   REAL(dp)                                :: dn_ppdzmin = 6._dp 
     274   REAL(dp)                                :: dn_pphmax  = 5750._dp 
     275   INTEGER(i4)                             :: in_nlevel  = 75 
    273276 
    274277   !namzps 
    275    REAL(dp)          :: dn_e3zps_min = 25._dp 
    276    REAL(dp)          :: dn_e3zps_rat = 0.2_dp 
     278   REAL(dp)                                :: dn_e3zps_min = 25._dp 
     279   REAL(dp)                                :: dn_e3zps_rat = 0.2_dp 
    277280 
    278281   ! namvar 
     282   CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varfile = '' 
    279283   CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varinfo = '' 
    280    CHARACTER(LEN=lc), DIMENSION(ip_maxvar) :: cn_varfile = '' 
    281284 
    282285   ! namnst 
    283    INTEGER(i4)       :: in_rhoi = 0 
    284    INTEGER(i4)       :: in_rhoj = 0 
     286   INTEGER(i4)                             :: in_rhoi = 0 
     287   INTEGER(i4)                             :: in_rhoj = 0 
    285288 
    286289   ! namout 
    287    CHARACTER(LEN=lc) :: cn_fileout = 'restart.nc'  
    288    LOGICAL           :: ln_extrap  = .FALSE. 
    289    INTEGER(i4)       :: in_nproc   = 0 
    290    INTEGER(i4)       :: in_niproc  = 0 
    291    INTEGER(i4)       :: in_njproc  = 0 
    292    CHARACTER(LEN=lc) :: cn_type    = '' 
     290   CHARACTER(LEN=lc)                       :: cn_fileout = 'restart.nc'  
     291   LOGICAL                                 :: ln_extrap  = .FALSE. 
     292   INTEGER(i4)                             :: in_nproc   = 0 
     293   INTEGER(i4)                             :: in_niproc  = 0 
     294   INTEGER(i4)                             :: in_njproc  = 0 
     295   CHARACTER(LEN=lc)                       :: cn_type    = '' 
    293296 
    294297   !------------------------------------------------------------------- 
     
    300303 
    301304   NAMELIST /namcfg/ &  !< configuration namelist 
    302    &  cn_varcfg         !< variable configuration file 
     305   &  cn_varcfg, &      !< variable configuration file 
     306   &  cn_dumcfg         !< dummy configuration file 
    303307 
    304308   NAMELIST /namcrs/ &  !< coarse grid namelist 
     
    330334 
    331335   NAMELIST /namvar/ &  !< variable namelist 
    332    &  cn_varinfo, &     !< list of variable and interpolation method to be used. 
    333    &  cn_varfile        !< list of variable file 
     336   &  cn_varfile, &     !< list of variable file 
     337   &  cn_varinfo        !< list of variable and interpolation method to be used. 
    334338    
    335339   NAMELIST /namnst/ &  !< nesting namelist 
     
    382386      ! get variable extra information 
    383387      CALL var_def_extra(TRIM(cn_varcfg)) 
     388 
     389      ! get dummy variable 
     390      CALL var_get_dummy(TRIM(cn_dumcfg)) 
     391      ! get dummy dimension 
     392      CALL dim_get_dummy(TRIM(cn_dumcfg)) 
     393      ! get dummy attribute 
     394      CALL att_get_dummy(TRIM(cn_dumcfg)) 
    384395 
    385396      READ( il_fileid, NML = namcrs ) 
     
    509520 
    510521               jvar=jvar+1 
    511                 
     522 
    512523               WRITE(*,'(2x,a,a)') "work on variable "//& 
    513524               &  TRIM(tl_multi%t_mpp(ji)%t_proc(1)%t_var(jj)%c_name) 
     
    541552            CALL iom_mpp_open(tl_mpp) 
    542553 
    543  
    544554            ! get or check depth value 
    545555            CALL create_restart_check_depth( tl_mpp, tl_depth ) 
     
    551561            CALL iom_mpp_close(tl_mpp) 
    552562 
    553             IF( ANY( tl_mpp%t_dim(1:2)%i_len /= & 
    554             &        tl_coord0%t_dim(1:2)%i_len) )THEN 
     563            IF( ANY(tl_mpp%t_dim(1:2)%i_len /= tl_coord0%t_dim(1:2)%i_len) .OR.& 
     564            &   ALL(il_rho(:)==1) )THEN 
    555565            !!! extract value from fine grid  
    556566 
    557                IF( ANY( tl_mpp%t_dim(1:2)%i_len <= & 
     567               IF( ANY( tl_mpp%t_dim(1:2)%i_len < & 
    558568               &        tl_coord1%t_dim(1:2)%i_len) )THEN 
    559                   CALL logger_fatal("CREATE RESTART: dimension in file "//& 
     569                  CALL logger_fatal("CREATE RESTART: dimensions in file "//& 
    560570                  &  TRIM(tl_mpp%c_name)//" smaller than those in fine"//& 
    561571                  &  " grid coordinates.") 
    562572               ENDIF 
    563573 
     574               ! use coord0 instead of mpp for restart file case  
     575               !  (without lon,lat) 
     576               ll_sameGrid=.FALSE. 
     577               IF( ALL(tl_mpp%t_dim(1:2)%i_len /= tl_coord0%t_dim(1:2)%i_len) & 
     578               &   )THEN 
     579                  ll_sameGrid=.TRUE.  
     580               ENDIF 
     581 
    564582               ! compute domain on fine grid 
    565                il_ind(:,:)=grid_get_coarse_index(tl_mpp, tl_coord1 ) 
     583               IF( ll_sameGrid )THEN 
     584                  il_ind(:,:)=grid_get_coarse_index(tl_mpp, tl_coord1 ) 
     585               ELSE 
     586                  il_ind(:,:)=grid_get_coarse_index(tl_coord0, tl_coord1 ) 
     587               ENDIF 
    566588 
    567589               il_imin1=il_ind(1,1) ; il_imax1=il_ind(1,2) 
     
    569591 
    570592               !- check grid coincidence 
    571                CALL grid_check_coincidence( tl_mpp, tl_coord1, & 
    572                &                            il_imin1, il_imax1, & 
    573                &                            il_jmin1, il_jmax1, & 
    574                &                            il_rho(:) ) 
     593               IF( ll_sameGrid )THEN 
     594                  CALL grid_check_coincidence( tl_mpp, tl_coord1, & 
     595                  &                            il_imin1, il_imax1, & 
     596                  &                            il_jmin1, il_jmax1, & 
     597                  &                            il_rho(:) ) 
     598               ELSE 
     599                  CALL grid_check_coincidence( tl_coord0, tl_coord1, & 
     600                  &                            il_imin1, il_imax1, & 
     601                  &                            il_jmin1, il_jmax1, & 
     602                  &                            il_rho(:) ) 
     603               ENDIF 
    575604 
    576605               ! compute domain 
     
    754783 
    755784   DO ji=1,ip_maxdim 
     785 
    756786      IF( tl_dim(ji)%l_use )THEN 
    757787         CALL mpp_move_dim(tl_mppout, tl_dim(ji)) 
     
    763793         END SELECT  
    764794      ENDIF 
     795 
    765796   ENDDO 
    766797 
     
    879910   !> and with dimension of the coordinate file.<br/>  
    880911   !> Then the variable array of value is split into equal subdomain. 
    881    !> Each subdomain is filled with the corresponding value of the matrix. 
     912   !> Each subdomain is filled with the associated value of the matrix. 
    882913   !> 
    883914   !> @author J.Paul 
     
    11691200            &        tl_depth%d_value(:,:,:,:) ) )THEN 
    11701201 
    1171                CALL logger_fatal("CREATE BOUNDARY: depth value from "//& 
    1172                &  TRIM(tl_multi%t_mpp(ji)%c_name)//" not conform "//& 
     1202               CALL logger_warn("CREATE BOUNDARY: depth value from "//& 
     1203               &  TRIM(td_mpp%c_name)//" not conform "//& 
    11731204               &  " to those from former file(s).") 
    11741205 
     
    12261257            IF( tl_date1 - tl_date2 /= 0 )THEN 
    12271258 
    1228                CALL logger_fatal("CREATE BOUNDARY: date from "//& 
    1229                &  TRIM(tl_multi%t_mpp(ji)%c_name)//" not conform "//& 
     1259               CALL logger_warn("CREATE BOUNDARY: date from "//& 
     1260               &  TRIM(td_mpp%c_name)//" not conform "//& 
    12301261               &  " to those from former file(s).") 
    12311262 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/dimension.f90

    r5617 r7339  
    154154! REVISION HISTORY: 
    155155!> @date November, 2013 - Initial Version 
     156!> @date Spetember, 2015 
     157!> - manage useless (dummy) dimension 
    156158!> 
    157159!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    167169   ! type and variable 
    168170   PUBLIC :: TDIM              !< dimension structure 
     171 
     172   PRIVATE :: cm_dumdim        !< dummy dimension array 
    169173 
    170174   ! function and subroutine 
     
    182186   PUBLIC :: dim_get_index     !< get dimension index in array of dimension structure 
    183187   PUBLIC :: dim_get_id        !< get dimension id in array of dimension structure 
     188   PUBLIC :: dim_get_dummy     !< fill dummy dimension array 
     189   PUBLIC :: dim_is_dummy      !< check if dimension is defined as dummy dimension 
    184190 
    185191   PRIVATE :: dim__reshape_2xyzt_dp ! reshape real(8) 4D array to ('x','y','z','t') 
     
    209215   END TYPE 
    210216 
     217   CHARACTER(LEN=lc), DIMENSION(ip_maxdum), SAVE :: cm_dumdim !< dummy dimension 
     218 
    211219   INTERFACE dim_print 
    212220      MODULE PROCEDURE dim__print_unit ! print information on one dimension 
     
    518526   !> @param[in] ld_uld    dimension unlimited 
    519527   !> @param[in] cd_sname  dimension short name 
    520    !> @param[in] ld_uld    dimension use or not 
     528   !> @param[in] ld_use    dimension use or not 
    521529   !> @return dimension structure 
    522530   !------------------------------------------------------------------- 
    523    TYPE(TDIM) FUNCTION dim_init( cd_name, id_len, ld_uld, cd_sname, ld_use) 
     531   TYPE(TDIM) FUNCTION dim_init( cd_name, id_len, ld_uld, cd_sname, ld_use ) 
    524532      IMPLICIT NONE 
    525533 
     
    14011409 
    14021410   END SUBROUTINE dim__clean_arr 
     1411   !------------------------------------------------------------------- 
     1412   !> @brief This subroutine fill dummy dimension array 
     1413   ! 
     1414   !> @author J.Paul 
     1415   !> @date September, 2015 - Initial Version 
     1416   ! 
     1417   !> @param[in] cd_dummy dummy configuration file 
     1418   !------------------------------------------------------------------- 
     1419   SUBROUTINE dim_get_dummy( cd_dummy ) 
     1420      IMPLICIT NONE 
     1421      ! Argument 
     1422      CHARACTER(LEN=*), INTENT(IN) :: cd_dummy 
     1423 
     1424      ! local variable 
     1425      INTEGER(i4)   :: il_fileid 
     1426      INTEGER(i4)   :: il_status 
     1427 
     1428      LOGICAL       :: ll_exist 
     1429 
     1430      ! loop indices 
     1431      ! namelist 
     1432      CHARACTER(LEN=lc), DIMENSION(ip_maxdum) :: cn_dumvar 
     1433      CHARACTER(LEN=lc), DIMENSION(ip_maxdum) :: cn_dumdim 
     1434      CHARACTER(LEN=lc), DIMENSION(ip_maxdum) :: cn_dumatt 
     1435 
     1436      !---------------------------------------------------------------- 
     1437      NAMELIST /namdum/ &   !< dummy namelist 
     1438      &  cn_dumvar, &       !< variable  name 
     1439      &  cn_dumdim, &       !< dimension name 
     1440      &  cn_dumatt          !< attribute name 
     1441      !---------------------------------------------------------------- 
     1442 
     1443      ! init 
     1444      cm_dumdim(:)='' 
     1445 
     1446      ! read namelist 
     1447      INQUIRE(FILE=TRIM(cd_dummy), EXIST=ll_exist) 
     1448      IF( ll_exist )THEN 
     1449 
     1450         il_fileid=fct_getunit() 
     1451 
     1452         OPEN( il_fileid, FILE=TRIM(cd_dummy), & 
     1453         &                FORM='FORMATTED',       & 
     1454         &                ACCESS='SEQUENTIAL',    & 
     1455         &                STATUS='OLD',           & 
     1456         &                ACTION='READ',          & 
     1457         &                IOSTAT=il_status) 
     1458         CALL fct_err(il_status) 
     1459         IF( il_status /= 0 )THEN 
     1460            CALL logger_fatal("DIM GET DUMMY: opening "//TRIM(cd_dummy)) 
     1461         ENDIF 
     1462 
     1463         READ( il_fileid, NML = namdum ) 
     1464         cm_dumdim(:)=cn_dumdim(:) 
     1465 
     1466         CLOSE( il_fileid ) 
     1467 
     1468      ENDIF 
     1469 
     1470   END SUBROUTINE dim_get_dummy 
     1471   !------------------------------------------------------------------- 
     1472   !> @brief This function check if dimension is defined as dummy dimension 
     1473   !> in configuraton file 
     1474   !> 
     1475   !> @author J.Paul 
     1476   !> @date September, 2015 - Initial Version 
     1477   ! 
     1478   !> @param[in] td_dim dimension structure 
     1479   !> @return true if dimension is dummy dimension  
     1480   !------------------------------------------------------------------- 
     1481   FUNCTION dim_is_dummy(td_dim) 
     1482      IMPLICIT NONE 
     1483 
     1484      ! Argument       
     1485      TYPE(TDIM), INTENT(IN) :: td_dim 
     1486       
     1487      ! function 
     1488      LOGICAL :: dim_is_dummy 
     1489       
     1490      ! loop indices 
     1491      INTEGER(i4) :: ji 
     1492      !---------------------------------------------------------------- 
     1493 
     1494      dim_is_dummy=.FALSE. 
     1495      DO ji=1,ip_maxdum 
     1496         IF( fct_lower(td_dim%c_name) == fct_lower(cm_dumdim(ji)) )THEN 
     1497            dim_is_dummy=.TRUE. 
     1498            EXIT 
     1499         ENDIF 
     1500      ENDDO 
     1501 
     1502   END FUNCTION dim_is_dummy 
    14031503END MODULE dim 
    14041504 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/docsrc/1_install.md

    r5617 r7339  
    1 # How to Install 
     1# Download 
    22 
    3 # Install NEMO 
    4 to install SIREN, you should first install NEMO. 
    5 see [here](http://www.nemo-ocean.eu/Using-NEMO/User-Guides/Basics/NEMO-Quick-Start-Guide) 
     3# Download NEMO # 
     4to install SIREN, you should first download NEMO. 
     5see [NEMO quick start guide](http://www.nemo-ocean.eu/Using-NEMO/User-Guides/Basics/NEMO-Quick-Start-Guide) 
    66 
    7 # Compile SIREN 
     7# Compile SIREN # 
    88when NEMO is installed, you just have to compile SIREN codes: 
    9 1. go to ./NEMOGCM/TOOLS 
    10 2. use maketools <br/> 
    11    to get help: maketools -h  
     9   1. go to ./NEMOGCM/TOOLS 
     10   2. run maketools (ex: ./maketools -n SIREN -m ifort_mpi_beaufix) 
    1211 
    13 # Fortran Compiler 
    14    SIREN codes were succesfully tested with : 
    15    - ifort (version 15.0.1) 
    16    - gfortran (version 4.8.2 20140120)  
    17 <!--   - pgf95 (version 13.9-0) --> 
     12      @note to get help on maketools: ./maketools -h 
    1813 
    19  <HR> 
    20    <b> 
    21    - @ref index 
    22    - @ref md_docsrc_3_codingRules 
    23    - @ref md_docsrc_4_changeLog 
    24    - @ref todo 
    25    </b> 
     14# Fortran Compiler # 
     15SIREN codes were succesfully tested with : 
     16  - ifort (version 15.0.1) 
     17  - gfortran (version 4.8.2 20140120)  
     18 
     19<HR> 
     20  <b> 
     21  - @ref index 
     22  - @ref md_docsrc_2_quickstart 
     23  - @ref md_docsrc_3_support_bug 
     24  - @ref md_docsrc_4_codingRules 
     25  - @ref md_docsrc_5_changeLog 
     26  - @ref todo 
     27  </b> 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/docsrc/main.dox

    r5037 r7339  
    11/*! 
    2  @mainpage Main Page 
    3  @section descr Generic Description 
    4  SIREN is a software to create regional configuration with 
    5  [NEMO](http://www.nemo-ocean.eu).<br/>  
     2 @mainpage About  
     3 
     4 SIREN is a software to create regional configuration with [NEMO](http://www.nemo-ocean.eu).<br/>  
    65 Actually SIREN create input files needed for a basic NEMO configuration.<br/> 
     6 
     7 SIREN allows you to create your own regional configuration embedded in a wider one.<br/> 
     8 In order to help you, a set of GLORYS files (global reanalysis on ORCA025 grid), as well as examples 
     9 of namelists are available in dods repository. 
     10 
     11 @note This software was created, and is maintain by the Configuration Manager Working Group, composed 
     12 of NEMO system team members. 
    713  
    8  SIREN is composed of a set of 5 Fortran programs : 
    9    - create_coord.f90 to create fine grid coordinate file from coarse grid coordinate file. 
    10    - create_bathy.f90 to create fine grid bathymetry file over domain. 
    11    - merge_bathy.f90 to merge fine grid bathymetry with coarse grid bathymetry at boundaries. 
    12    - create_restart.f90 to create initial state file from coarse grid restart or standard outputs. 
    13    - create_boundary.f90 to create boundary condition from coarse grid standard outputs. 
     14 To know how to install SIREN see @ref md_docsrc_1_install. 
    1415 
    15 To install those programs see @ref md_docsrc_1_install. 
    16  
    17  @note SIREN can not: 
    18  - create global configuration 
    19  - create configuarion around or close to north pole 
    20  - change number of vertical level 
    21  - change grid (horizontal or vertical) 
    22  
    23  @section howto How to use 
    24    @subsection howto_coord to create fine grid coordinate file 
    25    see create_coord.f90 
    26    @subsection howto_bathy to create fine grid bathymetry 
    27    see create_bathy.f90 
    28    @subsection howto_merge to merge fine grid bathymetry 
    29    see merge_bathy.f90 
    30    @subsection howto_restart to create initial state file 
    31    see create_restart.f90 
    32    @subsection howto_boundary to create boundary condition 
    33    see create_boundary.f90 
     16 You could find a tutorial for a quick start with SIREN in @ref md_docsrc_2_quickstart.<br/> 
     17 For more information about how to use each component of SIREN 
     18 - see create_coord.f90 to create fine grid coordinate file 
     19 - see create_bathy.f90 to create fine grid bathymetry 
     20 - see merge_bathy.f90 to merge fine grid bathymetry 
     21 - see create_restart.f90 to create initial state file, or other fields. 
     22 - see create_boundary.F90 to create boundary condition 
    3423 
    3524<HR> 
    3625   <b> 
    3726   - @ref md_docsrc_1_install 
    38    - @ref md_docsrc_3_codingRules 
    39    - @ref md_docsrc_4_changeLog 
     27   - @ref md_docsrc_2_quickstart 
     28   - @ref md_docsrc_3_support_bug 
     29   - @ref md_docsrc_4_codingRules 
     30   - @ref md_docsrc_5_changeLog 
    4031   - @ref todo 
    4132   </b> 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/domain.f90

    r5617 r7339  
    12971297   !> @date September, 2014 
    12981298   !> - take into account number of ghost cell 
     1299   !> @date February, 2016 
     1300   !> - number of extra point is the MAX (not the MIN) of zero and asess value.  
    12991301   ! 
    13001302   !> @param[inout] td_dom domain strcuture 
     
    13441346                  td_dom%i_imin      = td_dom%i_imin - td_dom%i_iextra(1) 
    13451347               ELSE ! td_dom%i_imin - il_iext <= td_dom%i_ghost0(jp_I,1)*ip_ghost 
    1346                   td_dom%i_iextra(1) = MIN(0, & 
     1348                  td_dom%i_iextra(1) = MAX(0, & 
    13471349                  &                         td_dom%i_imin - & 
    13481350                  &                         td_dom%i_ghost0(jp_I,1)*ip_ghost -1) 
     
    13561358               ELSE ! td_dom%i_imax + il_iext >= & 
    13571359                    !  td_dom%t_dim0(1)%i_len - td_dom%i_ghost0(jp_I,2)*ip_ghost 
    1358                   td_dom%i_iextra(2) = MIN(0, & 
     1360                  td_dom%i_iextra(2) = MAX( 0, & 
    13591361                  &                         td_dom%t_dim0(1)%i_len - & 
    13601362                  &                         td_dom%i_ghost0(jp_I,2)*ip_ghost - & 
     
    13641366 
    13651367            ELSE ! td_dom%i_ew0 >= 0 
     1368 
    13661369               ! EW cyclic 
    13671370               IF( td_dom%i_imin - il_iext > 0 )THEN 
     
    13911394            ! nothing to be done 
    13921395         ELSE 
     1396 
    13931397            IF( td_dom%i_jmin - il_jext > td_dom%i_ghost0(jp_J,1)*ip_ghost )THEN 
    13941398               td_dom%i_jextra(1) = il_jext 
    13951399               td_dom%i_jmin      = td_dom%i_jmin - td_dom%i_jextra(1) 
    13961400            ELSE ! td_dom%i_jmin - il_jext <= td_dom%i_ghost0(jp_J,1)*ip_ghost 
    1397                td_dom%i_jextra(1) = MIN(0, & 
     1401               td_dom%i_jextra(1) = MAX( 0, & 
    13981402               &                         td_dom%i_jmin - & 
    13991403               &                         td_dom%i_ghost0(jp_J,1)*ip_ghost - 1) 
     
    14071411            ELSE ! td_dom%i_jmax + il_jext >= & 
    14081412                 !  td_dom%t_dim0(2)%i_len - td_dom%i_ghost0(jp_J,2)*ip_ghost 
    1409                td_dom%i_jextra(2) = MIN(0, & 
     1413               td_dom%i_jextra(2) = MAX( 0, & 
    14101414               &                         td_dom%t_dim0(2)%i_len - & 
    14111415               &                         td_dom%i_ghost0(jp_J,2)*ip_ghost - & 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/file.f90

    r5617 r7339  
    694694   !> @date November, 2013 - Initial Version 
    695695   !> @date September, 2014 
    696    !> - add dimension to file if need be 
     696   !> - add dimension in file if need be 
    697697   !> - do not reorder dimension from variable, before put in file 
     698   !> @date September, 2015 
     699   !> - check variable dimension expected 
    698700   ! 
    699701   !> @param[inout] td_file   file structure 
     
    705707      ! Argument       
    706708      TYPE(TFILE), INTENT(INOUT) :: td_file 
    707       TYPE(TVAR) , INTENT(IN   ) :: td_var 
     709      TYPE(TVAR) , INTENT(INOUT) :: td_var 
    708710 
    709711      ! local variable 
     
    761763               IF( file_check_var_dim(td_file, td_var) )THEN 
    762764 
     765                  ! check variable dimension expected 
     766                  CALL var_check_dim(td_var) 
     767 
    763768                  ! update dimension if need be 
    764769                  DO ji=1,ip_maxdim 
     
    10501055                  ! new number of variable in file 
    10511056                  td_file%i_nvar=td_file%i_nvar-1 
    1052  
    10531057                  SELECT CASE(td_var%i_ndim) 
    10541058                     CASE(0) 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/function.f90

    r5609 r7339  
    363363      IF( id_status /= 0 )THEN 
    364364         !CALL ERRSNS() ! not F95 standard 
    365          PRINT *, "FORTRAN ERROR" 
     365         PRINT *, "FORTRAN ERROR ",id_status 
    366366         !STOP 
    367367      ENDIF 
     
    740740   ! 
    741741   !> @param[in] cd_var character 
    742    !> @return character is numeric 
     742   !> @return character is real number 
    743743   !------------------------------------------------------------------- 
    744744   PURE LOGICAL FUNCTION fct_is_real(cd_var) 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/global.f90

    r5037 r7339  
    1212! REVISION HISTORY: 
    1313!> @date November, 2013 - Initial Version 
     14!> @date September, 2015 
     15!> - define fill value for each variable type 
    1416! 
    1517!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    9597   &     'gauss      '/) 
    9698 
    97    REAL(dp)                                , PARAMETER :: dp_fill=NF90_FILL_DOUBLE !< default fill value 
     99   REAL(dp)                                , PARAMETER :: dp_fill_i1=NF90_FILL_BYTE   !< byte fill value 
     100   REAL(dp)                                , PARAMETER :: dp_fill_i2=NF90_FILL_SHORT  !< short fill value 
     101   REAL(dp)                                , PARAMETER :: dp_fill_i4=NF90_FILL_INT    !< INT fill value 
     102   REAL(dp)                                , PARAMETER :: dp_fill_sp=NF90_FILL_FLOAT  !< real fill value 
     103   REAL(dp)                                , PARAMETER :: dp_fill=NF90_FILL_DOUBLE !< double fill value 
    98104 
    99105   INTEGER(i4)                             , PARAMETER :: ip_npoint=4 
     
    125131   INTEGER(i4), PARAMETER :: jp_west =4 
    126132 
    127  
     133   INTEGER(i4)                             , PARAMETER :: ip_maxdum = 10 !< maximum dummy variable, dimension, attribute 
    128134 
    129135END MODULE global 
  • branches/2016/dev_NOC_2016/NEMOGCM/TOOLS/SIREN/src/grid.f90

    r5617 r7339  
    8080!> point:<br/> 
    8181!> @code 
    82 !>    il_index(:)=grid_get_closest(dd_lon0(:,:), dd_lat0(:,:), dd_lon1, dd_lat1) 
     82!>    il_index(:)=grid_get_closest(dd_lon0(:,:), dd_lat0(:,:), dd_lon1, dd_lat1 
     83!>                                 [,dd_fill] [,cd_pos]) 
    8384!> @endcode 
    8485!>       - il_index(:) is  coarse grid indices (/ i0, j0 /) 
     
    8788!>       - dd_lon1 is fine grid longitude value (real(8)) 
    8889!>       - dd_lat1 is fine grid latitude  value (real(8)) 
     90!>       - dd_fill 
     91!>       - cd_pos 
    8992!> 
    9093!>    to compute distance between a point A and grid points:<br/> 
     
    215218!> @date February, 2015 
    216219!> - add function grid_fill_small_msk to fill small domain inside bigger one 
     220!> @February, 2016 
     221!> - improve way to check coincidence (bug fix) 
     222!> - manage grid cases for T,U,V or F point, with even or odd refinment (bug fix) 
    217223! 
    218224!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    664670             
    665671            ! no pivot point found 
    666             CALL logger_error("GRID GET PIVOT: something wrong "//& 
     672            CALL logger_warn("GRID GET PIVOT: something wrong "//& 
    667673            &  "when computing pivot point with variable "//& 
    668674            &  TRIM(td_var%c_name)) 
     
    685691 
    686692               IF( grid__get_pivot_var /= -1 )THEN 
    687                   CALL logger_warn("GRID GET PIVOT: variable "//& 
     693                  CALL logger_info("GRID GET PIVOT: variable "//& 
    688694                  &  TRIM(td_var%c_name)//" seems to be on grid point "//& 
    689695                  &  TRIM(cp_grid_point(jj)) ) 
     
    13351341         il_dim(:)=td_var%t_dim(:)%i_len 
    13361342 
    1337          CALL logger_info("GRID GET PERIO: use varibale "//TRIM(td_var%c_name)) 
    1338          CALL logger_info("GRID GET PERIO: fill value "//TRIM(fct_str(td_var%d_fill))) 
    1339          CALL logger_info("GRID GET PERIO: fill value "//TRIM(fct_str(td_var%d_value(1,1,1,1)))) 
     1343         CALL logger_debug("GRID GET PERIO: use varibale "//TRIM(td_var%c_name)) 
     1344         CALL logger_debug("GRID GET PERIO: fill value "//TRIM(fct_str(td_var%d_fill))) 
     1345         CALL logger_debug("GRID GET PERIO: first value "//TRIM(fct_str(td_var%d_value(1,1,1,1)))) 
    13401346 
    13411347         IF(ALL(td_var%d_value(    1    ,    :    ,1,1)/=td_var%d_fill).AND.& 
     
    13441350         &  ALL(td_var%d_value(    :    ,il_dim(2),1,1)/=td_var%d_fill))THEN 
    13451351         ! no boundary closed 
    1346             CALL logger_warn("GRID GET PERIO: can't determined periodicity. "//& 
     1352            CALL logger_error("GRID GET PERIO: can't determined periodicity. "//& 
    13471353            &             "there is no boundary closed for variable "//& 
    13481354            &              TRIM(td_var%c_name) ) 
     1355            ! check pivot 
     1356            SELECT CASE(id_pivot) 
     1357               CASE(0) 
     1358                  ! F pivot 
     1359                  CALL logger_warn("GRID GET PERIO: assume domain is global") 
     1360                  grid__get_perio_var=6 
     1361               CASE(1) 
     1362                  ! T pivot 
     1363                  CALL logger_warn("GRID GET PERIO: assume domain is global") 
     1364                  grid__get_perio_var=4 
     1365            END SELECT 
    13491366         ELSE 
    13501367            ! check periodicity 
     
    22872304         &                                                    il_rho(:), cl_point ) 
    22882305 
    2289           
    22902306         CALL var_clean(tl_lon1) 
    22912307         CALL var_clean(tl_lat1)          
     
    24632479   !> - check grid point 
    24642480   !> - take into account EW overlap 
     2481   !> @date February, 2016 
     2482   !> - use delta (lon or lat) 
     2483   !> - manage cases for T,U,V or F point, with even or odd refinment 
    24652484   !> 
    24662485   !> @param[in] td_lon0   coarse grid longitude 
     
    24902509 
    24912510      ! local variable 
    2492       REAL(dp)    :: dl_lon1_ll 
    2493       REAL(dp)    :: dl_lon1_ul 
    2494       REAL(dp)    :: dl_lon1_lr 
    2495       REAL(dp)    :: dl_lon1_ur 
    2496  
    2497       REAL(dp)    :: dl_lat1_ll 
    2498       REAL(dp)    :: dl_lat1_ul 
    2499       REAL(dp)    :: dl_lat1_lr 
    2500       REAL(dp)    :: dl_lat1_ur 
     2511      CHARACTER(LEN= 1)                      :: cl_point0 
     2512      CHARACTER(LEN= 1)                      :: cl_point1 
     2513 
     2514      LOGICAL    , DIMENSION(2)              :: ll_even 
     2515 
     2516      REAL(dp)                               :: dl_lon1 
     2517      REAL(dp)                               :: dl_dlon 
     2518      REAL(dp)                               :: dl_lat1 
     2519      REAL(dp)                               :: dl_dlat 
     2520 
     2521      INTEGER(i4)                            :: il_ew0  
     2522      INTEGER(i4)                            :: il_imin0 
     2523      INTEGER(i4)                            :: il_imax0 
     2524      INTEGER(i4)                            :: il_jmin0 
     2525      INTEGER(i4)                            :: il_jmax0 
     2526 
     2527      INTEGER(i4)                            :: il_ew1  
     2528      INTEGER(i4)                            :: il_imin1 
     2529      INTEGER(i4)                            :: il_imax1 
     2530      INTEGER(i4)                            :: il_jmin1 
     2531      INTEGER(i4)                            :: il_jmax1 
     2532 
     2533      INTEGER(i4)                            :: il_imin 
     2534      INTEGER(i4)                            :: il_imax 
     2535      INTEGER(i4)                            :: il_jmin 
     2536      INTEGER(i4)                            :: il_jmax       
    25012537 
    25022538      INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho 
    25032539 
    2504       INTEGER(i4), DIMENSION(2) :: il_ill 
    2505       INTEGER(i4), DIMENSION(2) :: il_ilr 
    2506       INTEGER(i4), DIMENSION(2) :: il_iul 
    2507       INTEGER(i4), DIMENSION(2) :: il_iur 
    2508  
    2509       INTEGER(i4) :: il_ew0  
    2510       INTEGER(i4) :: il_imin0 
    2511       INTEGER(i4) :: il_imax0 
    2512       INTEGER(i4) :: il_jmin0 
    2513       INTEGER(i4) :: il_jmax0 
    2514  
    2515       INTEGER(i4) :: il_ew1  
    2516       INTEGER(i4) :: il_imin1 
    2517       INTEGER(i4) :: il_imax1 
    2518       INTEGER(i4) :: il_jmin1 
    2519       INTEGER(i4) :: il_jmax1 
    2520  
    2521       INTEGER(i4) :: il_imin 
    2522       INTEGER(i4) :: il_imax 
    2523       INTEGER(i4) :: il_jmin 
    2524       INTEGER(i4) :: il_jmax       
    2525  
    2526       INTEGER(i4), DIMENSION(2,2) :: il_xghost0 
    2527       INTEGER(i4), DIMENSION(2,2) :: il_yghost0 
    2528       INTEGER(i4), DIMENSION(2,2) :: il_xghost1 
    2529       INTEGER(i4), DIMENSION(2,2) :: il_yghost1 
    2530  
    2531       TYPE(TVAR) :: tl_lon0 
    2532       TYPE(TVAR) :: tl_lat0 
    2533       TYPE(TVAR) :: tl_lon1 
    2534       TYPE(TVAR) :: tl_lat1 
    2535  
    2536       CHARACTER(LEN= 1) :: cl_point0 
    2537       CHARACTER(LEN= 1) :: cl_point1 
    2538        
     2540      INTEGER(i4), DIMENSION(2)              :: il_ill 
     2541      INTEGER(i4), DIMENSION(2)              :: il_ilr 
     2542      INTEGER(i4), DIMENSION(2)              :: il_iul 
     2543      INTEGER(i4), DIMENSION(2)              :: il_iur 
     2544 
     2545      INTEGER(i4), DIMENSION(2,2)            :: il_xghost0 
     2546      INTEGER(i4), DIMENSION(2,2)            :: il_yghost0 
     2547      INTEGER(i4), DIMENSION(2,2)            :: il_xghost1 
     2548      INTEGER(i4), DIMENSION(2,2)            :: il_yghost1 
     2549 
     2550      TYPE(TVAR)                             :: tl_lon0 
     2551      TYPE(TVAR)                             :: tl_lat0 
     2552      TYPE(TVAR)                             :: tl_lon1 
     2553      TYPE(TVAR)                             :: tl_lat1 
     2554 
    25392555      ! loop indices 
    2540       INTEGER(i4) :: ji 
    2541       INTEGER(i4) :: jj 
    25422556      !---------------------------------------------------------------- 
    25432557      ! init 
     
    25472561      il_rho(:)=1 
    25482562      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:) 
     2563 
     2564      ll_even(:)=(/ (MOD(id_rho(jp_I),2)==0), (MOD(id_rho(jp_J),2)==0) /) 
    25492565 
    25502566      cl_point0='T' 
     
    26452661            ! get indices for each corner 
    26462662            !1- search lower left corner indices 
    2647             dl_lon1_ll=tl_lon1%d_value( il_imin1, il_jmin1, 1, 1 ) 
    2648             dl_lat1_ll=tl_lat1%d_value( il_imin1, il_jmin1, 1, 1 ) 
    2649  
    2650             IF( dl_lon1_ll == tl_lon1%d_fill .OR. & 
    2651             &   dl_lat1_ll == tl_lat1%d_fill )THEN 
    2652                CALL logger_debug("GRID GET COARSE INDEX: lon "//& 
    2653                &  TRIM(fct_str(dl_lon1_ll))//" "//& 
    2654                &  TRIM(fct_str(tl_lon1%d_fill)) ) 
    2655                CALL logger_debug("GRID GET COARSE INDEX: lat "//& 
    2656                &  TRIM(fct_str(dl_lat1_ll))//" "//& 
    2657                &  TRIM(fct_str(tl_lat1%d_fill)) ) 
     2663            dl_lon1=tl_lon1%d_value( il_imin1, il_jmin1, 1, 1 ) 
     2664            dl_lat1=tl_lat1%d_value( il_imin1, il_jmin1, 1, 1 ) 
     2665 
     2666            IF( dl_lon1 == tl_lon1%d_fill .OR. & 
     2667            &   dl_lat1 == tl_lat1%d_fill )THEN 
    26582668               CALL logger_error("GRID GET COARSE INDEX: lower left corner "//& 
    26592669               &                 "point is FillValue. remove ghost cell "//& 
    26602670               &                 "before running grid_get_coarse_index.") 
    26612671            ENDIF 
     2672 
     2673            !!!!! i-direction !!!!! 
     2674            IF( ll_even(jp_I) )THEN 
     2675               ! even 
     2676               SELECT CASE(TRIM(cl_point1)) 
     2677                  CASE('F','U') 
     2678                     dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmin1,1,1) -   & 
     2679                        &       tl_lon1%d_value(il_imin1  ,il_jmin1,1,1) ) / & 
     2680                        &     2. 
     2681                  CASE DEFAULT 
     2682                     dl_dlon=0 
     2683               END SELECT 
     2684            ELSE 
     2685               ! odd 
     2686               dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmin1,1,1) -   & 
     2687                  &       tl_lon1%d_value(il_imin1  ,il_jmin1,1,1) ) / & 
     2688                  &     2. 
     2689            ENDIF 
     2690 
     2691            !!!!! j-direction !!!!! 
     2692            IF( ll_even(jp_J) )THEN 
     2693               ! even 
     2694               SELECT CASE(TRIM(cl_point1)) 
     2695                  CASE('F','V') 
     2696                     dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmin1+1,1,1) -   & 
     2697                        &       tl_lat1%d_value(il_imin1,il_jmin1  ,1,1) ) / & 
     2698                        &     2. 
     2699                  CASE DEFAULT 
     2700                     dl_dlat=0 
     2701               END SELECT 
     2702            ELSE 
     2703               ! odd 
     2704               dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmin1+1,1,1) -   & 
     2705                  &       tl_lat1%d_value(il_imin1,il_jmin1  ,1,1) ) / & 
     2706                  &     2. 
     2707            ENDIF 
     2708 
     2709            dl_lon1 = dl_lon1 + dl_dlon 
     2710            dl_lat1 = dl_lat1 + dl_dlat 
     2711 
    26622712            ! look for closest point on coarse grid 
    26632713            il_ill(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, & 
     
    26672717            &                                           il_jmin0:il_jmax0, & 
    26682718            &                                           1,1), & 
    2669             &                           dl_lon1_ll, dl_lat1_ll   ) 
    2670  
    2671             ! coarse grid point should be south west of fine grid domain 
    2672             ji = il_ill(1) 
    2673             jj = il_ill(2) 
    2674  
    2675             IF( ABS(tl_lon0%d_value(ji,jj,1,1)-dl_lon1_ll) > dp_delta )THEN 
    2676                IF(tl_lon0%d_value(ji,jj,1,1) > dl_lon1_ll )THEN 
    2677                   il_ill(1)=il_ill(1)-1 
    2678                   IF( il_ill(1) <= 0 )THEN 
    2679                      IF( tl_lon0%i_ew >= 0 )THEN 
    2680                         il_ill(1)=tl_lon0%t_dim(jp_I)%i_len-tl_lon0%i_ew 
    2681                      ELSE 
    2682                         CALL logger_error("GRID GET COARSE INDEX: error "//& 
    2683                         &                 "computing lower left corner "//& 
    2684                         &                 "index for longitude") 
    2685                      ENDIF 
    2686                   ENDIF 
    2687                ENDIF 
    2688             ENDIF 
    2689             IF( ABS(tl_lat0%d_value(ji,jj,1,1)-dl_lat1_ll) > dp_delta )THEN 
    2690                IF(tl_lat0%d_value(ji,jj,1,1) > dl_lat1_ll )THEN 
    2691                   il_ill(2)=il_ill(2)-1 
    2692                   IF( il_ill(2)-1 <= 0 )THEN 
    2693                      CALL logger_error("GRID GET COARSE INDEX: error "//& 
    2694                      &                 "computing lower left corner "//& 
    2695                      &                 "index for latitude") 
    2696                   ENDIF 
    2697                ENDIF 
    2698             ENDIF 
     2719            &                           dl_lon1, dl_lat1, 'll'   ) 
     2720 
    26992721 
    27002722            !2- search upper left corner indices 
    2701             dl_lon1_ul=tl_lon1%d_value( il_imin1, il_jmax1, 1, 1 ) 
    2702             dl_lat1_ul=tl_lat1%d_value( il_imin1, il_jmax1, 1, 1 ) 
    2703  
    2704             IF( dl_lon1_ul == tl_lon1%d_fill .OR. & 
    2705             &   dl_lat1_ul == tl_lat1%d_fill )THEN 
     2723            dl_lon1=tl_lon1%d_value( il_imin1, il_jmax1, 1, 1 ) 
     2724            dl_lat1=tl_lat1%d_value( il_imin1, il_jmax1, 1, 1 ) 
     2725 
     2726            IF( dl_lon1 == tl_lon1%d_fill .OR. & 
     2727            &   dl_lat1 == tl_lat1%d_fill )THEN 
    27062728               CALL logger_error("GRID GET COARSE INDEX: upper left corner "//& 
    27072729               &                 "point is FillValue. remove ghost cell "//& 
    27082730               &                 "running grid_get_coarse_index.") 
    27092731            ENDIF             
     2732 
     2733            !!!!! i-direction !!!!! 
     2734            IF( ll_even(jp_I) )THEN 
     2735               ! even 
     2736               SELECT CASE(TRIM(cl_point1)) 
     2737                  CASE('F','U') 
     2738                     dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmax1,1,1) -   & 
     2739                        &       tl_lon1%d_value(il_imin1  ,il_jmax1,1,1) ) / & 
     2740                        &     2. 
     2741                  CASE DEFAULT 
     2742                     dl_dlon=0 
     2743               END SELECT 
     2744            ELSE 
     2745               ! odd 
     2746               dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmax1,1,1) -   & 
     2747                  &       tl_lon1%d_value(il_imin1  ,il_jmax1,1,1) ) / & 
     2748                  &     2. 
     2749            ENDIF 
     2750 
     2751            !!!!! j-direction !!!!! 
     2752            IF( ll_even(jp_J) )THEN 
     2753               ! even 
     2754               SELECT CASE(TRIM(cl_point1)) 
     2755                  CASE('F','V') 
     2756                     dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmax1  ,1,1) -   & 
     2757                        &       tl_lat1%d_value(il_imin1,il_jmax1-1,1,1) ) / & 
     2758                        &     2. 
     2759                  CASE DEFAULT 
     2760                     dl_dlat=0 
     2761               END SELECT 
     2762            ELSE 
     2763               ! odd 
     2764               dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmax1  ,1,1) -   & 
     2765                  &       tl_lat1%d_value(il_imin1,il_jmax1-1,1,1) ) / & 
     2766                  &     2. 
     2767            ENDIF 
     2768 
     2769            dl_lon1 = dl_lon1 + dl_dlon 
     2770            dl_lat1 = dl_lat1 - dl_dlat 
     2771 
    27102772            ! look for closest point on coarse grid 
    27112773            il_iul(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, & 
     
    27152777            &                                           il_jmin0:il_jmax0, & 
    27162778            &                                           1,1), & 
    2717             &                           dl_lon1_ul, dl_lat1_ul   ) 
    2718  
    2719             ! coarse grid point should be north west of fine grid domain 
    2720             ji = il_iul(1) 
    2721             jj = il_iul(2) 
    2722             IF( ABS(tl_lon0%d_value(ji,jj,1,1)-dl_lon1_ul) > dp_delta )THEN 
    2723                IF(tl_lon0%d_value(ji,jj,1,1) > dl_lon1_ul )THEN 
    2724                   il_iul(1)=il_iul(1)-1 
    2725                   IF( il_iul(1) <= 0 )THEN 
    2726                      IF( tl_lon0%i_ew >= 0 )THEN 
    2727                         il_iul(1)=tl_lon0%t_dim(jp_I)%i_len-tl_lon0%i_ew 
    2728                      ELSE 
    2729                         CALL logger_error("GRID GET COARSE INDEX: error "//& 
    2730                         &                 "computing upper left corner "//& 
    2731                         &                 "index for longitude") 
    2732                      ENDIF 
    2733                   ENDIF 
    2734                ENDIF 
    2735             ENDIF 
    2736             IF( ABS(tl_lat0%d_value(ji,jj,1,1)-dl_lat1_ul) > dp_delta )THEN 
    2737                IF(tl_lat0%d_value(ji,jj,1,1) < dl_lat1_ul )THEN 
    2738                   il_iul(2)=il_iul(2)+1 
    2739                   IF( il_ill(2) > tl_lat0%t_dim(jp_J)%i_len )THEN 
    2740                      CALL logger_error("GRID GET COARSE INDEX: error "//& 
    2741                      &                 "computing upper left corner "//& 
    2742                      &                 "index for latitude") 
    2743                   ENDIF 
    2744                ENDIF 
    2745             ENDIF 
     2779            &                           dl_lon1, dl_lat1, 'ul' ) 
    27462780 
    27472781            !3- search lower right corner indices 
    2748             dl_lon1_lr=tl_lon1%d_value( il_imax1, il_jmin1, 1, 1 ) 
    2749             dl_lat1_lr=tl_lat1%d_value( il_imax1, il_jmin1, 1, 1 ) 
    2750  
    2751             IF( dl_lon1_lr == tl_lon1%d_fill .OR. & 
    2752             &   dl_lat1_lr == tl_lat1%d_fill )THEN 
     2782            dl_lon1=tl_lon1%d_value( il_imax1, il_jmin1, 1, 1 ) 
     2783            dl_lat1=tl_lat1%d_value( il_imax1, il_jmin1, 1, 1 ) 
     2784 
     2785            IF( dl_lon1 == tl_lon1%d_fill .OR. & 
     2786            &   dl_lat1 == tl_lat1%d_fill )THEN 
    27532787               CALL logger_error("GRID GET COARSE INDEX: lower right corner "//& 
    27542788               &                 "point is FillValue. remove ghost cell "//& 
    27552789               &                 "running grid_get_coarse_index.") 
    27562790            ENDIF             
     2791 
     2792            !!!!! i-direction !!!!! 
     2793            IF( ll_even(jp_I) )THEN 
     2794               ! even 
     2795               SELECT CASE(TRIM(cl_point1)) 
     2796                  CASE('F','U') 
     2797                     dl_dlon= ( tl_lon1%d_value(il_imax1  ,il_jmin1,1,1) -   & 
     2798                        &       tl_lon1%d_value(il_imax1-1,il_jmin1,1,1) ) / & 
     2799                        &     2. 
     2800                  CASE DEFAULT 
     2801                     dl_dlon=0 
     2802               END SELECT 
     2803            ELSE 
     2804               ! odd 
     2805               dl_dlon= ( tl_lon1%d_value(il_imax1  ,il_jmin1,1,1) -   & 
     2806                  &       tl_lon1%d_value(il_imax1-1,il_jmin1,1,1) ) / & 
     2807                  &     2. 
     2808            ENDIF 
     2809 
     2810            !!!!! j-direction !!!!! 
     2811            IF( ll_even(jp_J) )THEN 
     2812               ! even 
     2813               SELECT CASE(TRIM(cl_point1)) 
     2814                  CASE('F','V') 
     2815                     dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmin1+1,1,1) -   & 
     2816                        &       tl_lat1%d_value(il_imax1,il_jmin1  ,1,1) ) / & 
     2817                        &     2. 
     2818                  CASE DEFAULT 
     2819                     dl_dlat=0 
     2820               END SELECT 
     2821            ELSE 
     2822               ! odd 
     2823               dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmin1+1,1,1) -   & 
     2824                  &       tl_lat1%d_value(il_imax1,il_jmin1  ,1,1) ) / & 
     2825                  &     2. 
     2826            ENDIF 
     2827 
     2828            dl_lon1 = dl_lon1 - dl_dlon 
     2829            dl_lat1 = dl_lat1 + dl_dlat 
     2830 
    27572831            ! look for closest point on coarse grid 
    27582832            il_ilr(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, & 
     
    27622836            &                                           il_jmin0:il_jmax0, & 
    27632837            &                                           1,1), & 
    2764             &                           dl_lon1_lr, dl_lat1_lr   ) 
    2765  
    2766             ! coarse grid point should be south east of fine grid domain 
    2767             ji = il_ilr(1) 
    2768             jj = il_ilr(2) 
    2769             IF( ABS(tl_lon0%d_value(ji,jj,1,1)-dl_lon1_lr) > dp_delta )THEN 
    2770                IF( tl_lon0%d_value(ji,jj,1,1) < dl_lon1_lr )THEN 
    2771                   il_ilr(1)=il_ilr(1)+1 
    2772                   IF( il_ilr(1) > tl_lon0%t_dim(jp_I)%i_len )THEN 
    2773                      IF( tl_lon0%i_ew >= 0 )THEN 
    2774                         il_ilr(1)=tl_lon0%i_ew+1 
    2775                      ELSE 
    2776                         CALL logger_error("GRID GET COARSE INDEX: error "//& 
    2777                         &                 "computing lower right corner "//& 
    2778                         &                 "index for longitude") 
    2779                      ENDIF 
    2780                   ENDIF 
    2781                ENDIF 
    2782             ENDIF 
    2783             IF( ABS(tl_lat0%d_value(ji,jj,1,1)-dl_lat1_lr) > dp_delta )THEN 
    2784                IF( tl_lat0%d_value(ji,jj,1,1) > dl_lat1_lr )THEN 
    2785                   il_ilr(2)=il_ilr(2)-1 
    2786                   IF( il_ilr(2) <= 0 )THEN 
    2787                      CALL logger_error("GRID GET COARSE INDEX: error "//& 
    2788                      &                 "computing lower right corner "//& 
    2789                      &                 "index for latitude") 
    2790                   ENDIF 
    2791                ENDIF 
    2792             ENDIF 
     2838            &                           dl_lon1, dl_lat1, 'lr' ) 
    27932839 
    27942840            !4- search upper right corner indices 
    2795             dl_lon1_ur=tl_lon1%d_value( il_imax1, il_jmax1, 1, 1 ) 
    2796             dl_lat1_ur=tl_lat1%d_value( il_imax1, il_jmax1, 1, 1 ) 
    2797  
    2798             IF( dl_lon1_ur == tl_lon1%d_fill .OR. & 
    2799             &   dl_lat1_ur == tl_lat1%d_fill )THEN 
     2841            dl_lon1=tl_lon1%d_value( il_imax1, il_jmax1, 1, 1 ) 
     2842            dl_lat1=tl_lat1%d_value( il_imax1, il_jmax1, 1, 1 ) 
     2843 
     2844            IF( dl_lon1 == tl_lon1%d_fill .OR. & 
     2845            &   dl_lat1 == tl_lat1%d_fill )THEN 
    28002846               CALL logger_error("GRID GET COARSE INDEX: upper right corner "//& 
    28012847               &                 "point is FillValue. remove ghost cell "//& 
    2802                &                 "running grid_get_coarse_index.") 
     2848               &                 "before running grid_get_coarse_index.") 
    28032849            ENDIF             
     2850 
     2851            !!!!! i-direction !!!!! 
     2852            IF( ll_even(jp_I) )THEN 
     2853               ! even 
     2854               SELECT CASE(TRIM(cl_point1)) 
     2855                  CASE('F','U') 
     2856                     dl_dlon= ( tl_lon1%d_value(il_imax1  ,il_jmax1,1,1) -   & 
     2857                        &       tl_lon1%d_value(il_imax1-1,il_jmax1,1,1) ) / & 
     2858                        &     2. 
     2859                  CASE DEFAULT 
     2860                     dl_dlon=0 
     2861               END SELECT 
     2862            ELSE 
     2863               ! odd 
     2864               dl_dlon= ( tl_lon1%d_value(il_imax1  ,il_jmax1,1,1) -   & 
     2865                  &       tl_lon1%d_value(il_imax1-1,il_jmax1,1,1) ) / & 
     2866                  &     2. 
     2867            ENDIF 
     2868 
     2869            !!!!! j-direction !!!!! 
     2870            IF( ll_even(jp_J) )THEN 
     2871               ! even 
     2872               SELECT CASE(TRIM(cl_point1)) 
     2873                  CASE('F','V') 
     2874                     dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmax1  ,1,1) -   & 
     2875                        &       tl_lat1%d_value(il_imax1,il_jmax1-1,1,1) ) / & 
     2876                        &     2. 
     2877                  CASE DEFAULT 
     2878                     dl_dlat=0 
     2879               END SELECT 
     2880            ELSE 
     2881               ! odd 
     2882               dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmax1  ,1,1) -   & 
     2883                  &       tl_lat1%d_value(il_imax1,il_jmax1-1,1,1) ) / & 
     2884                  &     2. 
     2885            ENDIF 
     2886 
     2887            dl_lon1 = dl_lon1 - dl_dlon 
     2888            dl_lat1 = dl_lat1 - dl_dlat 
     2889 
    28042890            ! look for closest point on coarse grid 
    28052891            il_iur(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, & 
     
    28092895            &                                           il_jmin0:il_jmax0, & 
    28102896            &                                           1,1), & 
    2811             &                           dl_lon1_ur, dl_lat1_ur   ) 
    2812  
    2813             ! coarse grid point should be north east fine grid domain 
    2814             ji = il_iur(1) 
    2815             jj = il_iur(2) 
    2816             IF( ABS(tl_lon0%d_value(ji,jj,1,1)-dl_lon1_ur) > dp_delta )THEN 
    2817                IF( tl_lon0%d_value(ji,jj,1,1) < dl_lon1_ur )THEN 
    2818                   il_iur(1)=il_iur(1)+1 
    2819                   IF( il_iur(1) > tl_lon0%t_dim(jp_I)%i_len )THEN 
    2820                      IF( tl_lon0%i_ew >= 0 )THEN 
    2821                         il_iur(1)=tl_lon0%i_ew+1 
    2822                      ELSE 
    2823                         CALL logger_error("GRID GET COARSE INDEX: error "//& 
    2824                         &                 "computing upper right corner "//& 
    2825                         &                 "index for longitude") 
    2826                      ENDIF 
    2827                   ENDIF 
    2828                ENDIF 
    2829             ENDIF 
    2830             IF( ABS(tl_lat0%d_value(ji,jj,1,1)-dl_lat1_ur) > dp_delta )THEN 
    2831                IF( tl_lat0%d_value(ji,jj,1,1) < dl_lat1_ur )THEN 
    2832                   il_iur(2)=il_iur(2)+1 
    2833                   IF( il_iur(2) > tl_lat0%t_dim(jp_J)%i_len )THEN 
    2834                      CALL logger_error("GRID GET COARSE INDEX: error "//& 
    2835                      &                 "computing upper right corner "//& 
    2836                      &                 "index for latitude") 
    2837                   ENDIF 
    2838                ENDIF 
    2839             ENDIF 
     2897            &                           dl_lon1, dl_lat1, 'ur' ) 
    28402898 
    28412899            ! coarse grid indices 
     
    29433001   END FUNCTION grid_is_global 
    29443002   !------------------------------------------------------------------- 
    2945    !> @brief This function return coarse grid indices of the closest point 
    2946    !> from fine grid point (lon1,lat1)  
     3003   !> @brief This function return grid indices of the closest point 
     3004   !> from point (lon1,lat1)  
    29473005   !>  
    29483006   !> @details 
     
    29513009   !> of longitude and latitude, before running this function 
    29523010   !> 
     3011   !> if you add cd_pos argument, you could choice to return closest point at 
     3012   !> - lower left  (ll) of the point 
     3013   !> - lower right (lr) of the point 
     3014   !> - upper left  (ul) of the point 
     3015   !> - upper right (ur) of the point 
     3016   !> - lower       (lo) of the point 
     3017   !> - upper       (up) of the point 
     3018   !> -       left  (le) of the point 
     3019   !> -       right (ri) of the point 
     3020   !> 
    29533021   !> @author J.Paul 
    29543022   !> @date November, 2013 - Initial Version 
    2955    !> @date February, 2015 - change dichotomy method to manage ORCA grid 
     3023   !> @date February, 2015 
     3024   !> - change dichotomy method to manage ORCA grid 
     3025   !> @date February, 2016 
     3026   !> - add optional use of relative position 
    29563027   ! 
    29573028   !> @param[in] dd_lon0   coarse grid array of longitude 
     
    29593030   !> @param[in] dd_lon1   fine   grid longitude 
    29603031   !> @param[in] dd_lat1   fine   grid latitude 
     3032   !> @param[in] cd_pos    relative position of grid point from point  
    29613033   !> @param[in] dd_fill   fill value 
    29623034   !> @return coarse grid indices of closest point of fine grid point 
    29633035   !------------------------------------------------------------------- 
    2964    FUNCTION grid_get_closest( dd_lon0, dd_lat0, dd_lon1, dd_lat1, dd_fill ) 
     3036   FUNCTION grid_get_closest( dd_lon0, dd_lat0, dd_lon1, dd_lat1, cd_pos, dd_fill ) 
    29653037      IMPLICIT NONE 
    29663038      ! Argument 
     
    29693041      REAL(dp),                 INTENT(IN) :: dd_lon1 
    29703042      REAL(dp),                 INTENT(IN) :: dd_lat1 
     3043      CHARACTER(LEN=*),         INTENT(IN), OPTIONAL :: cd_pos 
    29713044      REAL(dp),                 INTENT(IN), OPTIONAL :: dd_fill 
    29723045 
     
    31473220      &                          dl_lon1, dd_lat1 ) 
    31483221 
     3222      IF( PRESENT(cd_pos) )THEN 
     3223         !  
     3224         SELECT CASE(TRIM(cd_pos)) 
     3225            CASE('le') 
     3226               WHERE( dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) > dd_lat1 ) 
     3227                  dl_dist(:,:)=NF90_FILL_DOUBLE 
     3228               END WHERE 
     3229            CASE('ri') 
     3230               WHERE( dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) < dd_lat1 ) 
     3231                  dl_dist(:,:)=NF90_FILL_DOUBLE 
     3232               END WHERE 
     3233            CASE('up') 
     3234               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) < dl_lon1 ) 
     3235                  dl_dist(:,:)=NF90_FILL_DOUBLE 
     3236               END WHERE 
     3237            CASE('lo') 
     3238               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) > dl_lon1 ) 
     3239                  dl_dist(:,:)=NF90_FILL_DOUBLE 
     3240               END WHERE 
     3241            CASE('ll') 
     3242               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) > dl_lon1 .OR. & 
     3243                    & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) > dd_lat1 ) 
     3244                  dl_dist(:,:)=NF90_FILL_DOUBLE 
     3245               END WHERE 
     3246            CASE('lr') 
     3247               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) < dl_lon1 .OR. & 
     3248                    & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) > dd_lat1 ) 
     3249                  dl_dist(:,:)=NF90_FILL_DOUBLE 
     3250               END WHERE                
     3251            CASE('ul') 
     3252               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) > dl_lon1 .OR. & 
     3253                    & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) < dd_lat1 ) 
     3254                  dl_dist(:,:)=NF90_FILL_DOUBLE 
     3255               END WHERE                
     3256            CASE('ur') 
     3257               WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) < dl_lon1 .OR. & 
     3258                    & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) < dd_lat1 ) 
     3259                  dl_dist(:,:)=NF90_FILL_DOUBLE 
     3260               END WHERE 
     3261         END SELECT 
     3262      ENDIF 
    31493263      grid_get_closest(:)=MINLOC(dl_dist(:,:),dl_dist(:,:)/=NF90_FILL_DOUBLE) 
    31503264 
     
    34433557         &                                         il_imax0, il_jmax0, & 
    34443558         &                                         dl_lon1(:,:), dl_lat1(:,:),& 
    3445          &                                         id_rho(:) ) 
     3559         &                                         id_rho(:), cl_point ) 
    34463560  
    34473561         DEALLOCATE(dl_lon0, dl_lat0) 
     
    35883702         &                                         id_imax0, id_jmax0, & 
    35893703         &                                         dl_lon1(:,:), dl_lat1(:,:),& 
    3590          &                                         id_rho(:) ) 
     3704         &                                         id_rho(:), cl_point ) 
    35913705          
    35923706         DEALLOCATE(dl_lon1, dl_lat1) 
     
    36683782      ! init 
    36693783      grid__get_fine_offset_fc(:,:)=-1 
    3670  
    36713784      ALLOCATE(il_rho(ip_maxdim)) 
    36723785      il_rho(:)=1 
     
    36903803         CALL iom_mpp_open(tl_coord0) 
    36913804 
    3692          ! read coarse longitue and latitude 
     3805         ! read coarse longitude and latitude 
    36933806         WRITE(cl_name,*) 'longitude_'//TRIM(cl_point) 
    36943807         il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name) 
     
    37103823         ENDIF 
    37113824         tl_lat0=iom_mpp_read_var(tl_coord0, TRIM(cl_name)) 
    3712           
     3825  
    37133826         ! close mpp files 
    37143827         CALL iom_mpp_close(tl_coord0) 
     
    37163829         CALL grid_del_ghost(tl_lon0, il_xghost0(:,:)) 
    37173830         CALL grid_del_ghost(tl_lat0, il_xghost0(:,:)) 
     3831 
    37183832 
    37193833         ALLOCATE(dl_lon0(tl_lon0%t_dim(jp_I)%i_len, & 
     
    37383852         il_jmax0=id_jmax0-il_xghost0(jp_J,1) 
    37393853 
    3740        
    37413854         !3- compute 
    37423855         grid__get_fine_offset_fc(:,:)=grid_get_fine_offset(& 
     
    37453858         &                                         il_imax0, il_jmax0, & 
    37463859         &                                         dd_lon1(:,:), dd_lat1(:,:),& 
    3747          &                                         id_rho(:) ) 
     3860         &                                         id_rho(:), cl_point ) 
    37483861          
    37493862         DEALLOCATE(dl_lon0, dl_lat0) 
     
    37673880   !> @date May, 2015  
    37683881   !> - improve way to find offset 
     3882   !> @date July, 2015 
     3883   !> - manage case close to greenwich meridian 
     3884   !> @date February, 2016 
     3885   !> - use grid_get_closest to assess offset 
     3886   !> - use delta (lon or lat) 
     3887   !> - manage cases for T,U,V or F point, with even or odd refinment 
     3888   !> - check lower left(upper right) fine grid point inside lower left(upper 
     3889   !> right) coarse grid cell. 
     3890   !>  
     3891   !> @todo check case close from North fold. 
    37693892   !> 
    37703893   !> @param[in] dd_lon0   coarse grid longitude array  
     
    37773900   !> @param[in] dd_lat1   fine   grid latitude  array 
    37783901   !> @param[in] id_rho    array of refinement factor 
     3902   !> @param[in] cd_point  Arakawa grid point 
    37793903   !> @return offset array (/ (/i_offset_left,i_offset_right/),(/j_offset_lower,j_offset_upper/) /) 
    37803904   !------------------------------------------------------------------- 
    37813905   FUNCTION grid__get_fine_offset_cc( dd_lon0, dd_lat0, & 
    37823906   &                                  id_imin0, id_jmin0, id_imax0, id_jmax0, & 
    3783    &                                  dd_lon1, dd_lat1, id_rho ) 
     3907   &                                  dd_lon1, dd_lat1, id_rho, cd_point ) 
    37843908      IMPLICIT NONE 
    37853909      ! Argument 
    3786       REAL(dp)   , DIMENSION(:,:), INTENT(IN) :: dd_lon0 
    3787       REAL(dp)   , DIMENSION(:,:), INTENT(IN) :: dd_lat0 
    3788       REAL(dp)   , DIMENSION(:,:), INTENT(IN) :: dd_lon1 
    3789       REAL(dp)   , DIMENSION(:,:), INTENT(IN) :: dd_lat1 
    3790  
    3791       INTEGER(i4),                 INTENT(IN) :: id_imin0 
    3792       INTEGER(i4),                 INTENT(IN) :: id_jmin0 
    3793       INTEGER(i4),                 INTENT(IN) :: id_imax0 
    3794       INTEGER(i4),                 INTENT(IN) :: id_jmax0 
    3795  
    3796       INTEGER(i4), DIMENSION(:)  , INTENT(IN) :: id_rho 
     3910      REAL(dp)        , DIMENSION(:,:), INTENT(IN) :: dd_lon0 
     3911      REAL(dp)        , DIMENSION(:,:), INTENT(IN) :: dd_lat0 
     3912      REAL(dp)        , DIMENSION(:,:), INTENT(IN) :: dd_lon1 
     3913      REAL(dp)        , DIMENSION(:,:), INTENT(IN) :: dd_lat1 
     3914 
     3915      INTEGER(i4)     ,                 INTENT(IN) :: id_imin0 
     3916      INTEGER(i4)     ,                 INTENT(IN) :: id_jmin0 
     3917      INTEGER(i4)     ,                 INTENT(IN) :: id_imax0 
     3918      INTEGER(i4)     ,                 INTENT(IN) :: id_jmax0 
     3919 
     3920      INTEGER(i4)     , DIMENSION(:)  , INTENT(IN) :: id_rho 
     3921      CHARACTER(LEN=*)                , INTENT(IN), OPTIONAL :: cd_point 
    37973922 
    37983923      ! function 
     
    38003925 
    38013926      ! local variable 
     3927      CHARACTER(LEN= 1)                        :: cl_point 
     3928 
     3929      INTEGER(i4)                              :: i1 
     3930      INTEGER(i4)                              :: i2 
     3931      INTEGER(i4)                              :: j1 
     3932      INTEGER(i4)                              :: j2 
     3933 
    38023934      INTEGER(i4), DIMENSION(2)                :: il_shape0 
    38033935      INTEGER(i4), DIMENSION(2)                :: il_shape1 
    38043936 
     3937      INTEGER(i4), DIMENSION(2)                :: il_ind 
     3938 
    38053939      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon0 
    38063940      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_lon1 
    38073941 
    3808       LOGICAL                                  :: ll_ii 
    3809       LOGICAL                                  :: ll_ij 
     3942      REAL(dp)                                 :: dl_lonmax0 
     3943      REAL(dp)                                 :: dl_latmax0 
     3944      REAL(dp)                                 :: dl_lonmin0 
     3945      REAL(dp)                                 :: dl_latmin0 
     3946 
     3947      REAL(dp)                                 :: dl_lon0F 
     3948      REAL(dp)                                 :: dl_lat0F 
     3949      REAL(dp)                                 :: dl_dlon 
     3950      REAL(dp)                                 :: dl_dlat 
     3951 
     3952      LOGICAL    , DIMENSION(2)                :: ll_even 
     3953      LOGICAL                                  :: ll_greenwich 
    38103954       
    38113955      ! loop indices 
    3812       INTEGER(i4) :: ji 
    3813       INTEGER(i4) :: jj 
    3814  
    38153956      INTEGER(i4) :: ii 
    38163957      INTEGER(i4) :: ij 
     
    38243965         CALL logger_fatal("GRID GET FINE OFFSET: dimension of fine "//& 
    38253966         &              "longitude and latitude differ") 
    3826       ENDIF       
     3967      ENDIF 
     3968 
     3969      ll_even(:)=(/ (MOD(id_rho(jp_I),2)==0), (MOD(id_rho(jp_J),2)==0) /) 
     3970 
     3971      cl_point='T' 
     3972      IF( PRESENT(cd_point) ) cl_point=TRIM(fct_upper(cd_point)) 
    38273973 
    38283974      il_shape0(:)=SHAPE(dd_lon0(:,:)) 
    38293975      ALLOCATE( dl_lon0(il_shape0(1),il_shape0(2)) ) 
    38303976 
     3977      il_shape1(:)=SHAPE(dd_lon1(:,:)) 
     3978      ALLOCATE( dl_lon1(il_shape1(1),il_shape1(2)) ) 
     3979 
    38313980      dl_lon0(:,:)=dd_lon0(:,:) 
    38323981      WHERE( dd_lon0(:,:) < 0 ) dl_lon0(:,:)=dd_lon0(:,:)+360. 
    38333982 
    3834       il_shape1(:)=SHAPE(dd_lon1(:,:)) 
    3835       ALLOCATE( dl_lon1(il_shape1(1),il_shape1(2)) ) 
    3836  
    38373983      dl_lon1(:,:)=dd_lon1(:,:) 
    3838       WHERE( dd_lon1(:,:) < 0 ) dl_lon1(:,:)=dd_lon1(:,:)+360.          
     3984      WHERE( dd_lon1(:,:) < 0 ) dl_lon1(:,:)=dd_lon1(:,:)+360. 
    38393985 
    38403986      ! init 
    38413987      grid__get_fine_offset_cc(:,:)=-1 
     3988      ll_greenwich=.FALSE. 
    38423989 
    38433990      IF( il_shape1(jp_J) == 1 )THEN 
    3844            
     3991  
    38453992         grid__get_fine_offset_cc(jp_J,:)=((id_rho(jp_J)-1)/2) 
    38463993 
    3847          ! work on i-direction 
    3848          ! look for i-direction left offset 
    3849          IF( dl_lon1(1,1) < dl_lon0(id_imin0+1,id_jmin0) )THEN 
    3850             DO ji=1,id_rho(jp_I)+2 
    3851                IF( dl_lon1(ji,1) > dl_lon0(id_imin0+1,id_jmin0) - dp_delta )THEN 
    3852                   grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ji 
    3853                   EXIT 
    3854                ENDIF 
    3855             ENDDO 
     3994         !!! work on i-direction 
     3995         !!! look for i-direction left offset 
     3996         i1=1 ; i2=MIN((id_rho(jp_I)+2),il_shape1(jp_I)) 
     3997         j1=1 ; j2=1 
     3998 
     3999         ! check if cross greenwich meridien 
     4000         IF( minval(dl_lon0(id_imin0:id_imin0+1,id_jmin0))<5. .OR. & 
     4001           & maxval(dl_lon0(id_imin0:id_imin0+1,id_jmin0))>355. )THEN 
     4002            ! close to greenwich meridien 
     4003            ll_greenwich=.TRUE. 
     4004            ! 0:360 => -180:180 
     4005            WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0) > 180. ) 
     4006               dl_lon0(id_imin0:id_imin0+1,id_jmin0) = & 
     4007                  & dl_lon0(id_imin0:id_imin0+1,id_jmin0)-360. 
     4008            END WHERE 
     4009 
     4010            WHERE( dl_lon1(i1:i2,j1:j2) > 180. ) 
     4011               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360. 
     4012            END WHERE 
     4013         ENDIF 
     4014 
     4015         ! max lognitude of the left cell 
     4016         dl_lonmax0=dl_lon0(id_imin0+1,id_jmin0) 
     4017         IF( dl_lon1(1,1) < dl_lonmax0 )THEN 
     4018 
     4019            !!!!! i-direction !!!!! 
     4020            IF( ll_even(jp_I) )THEN 
     4021               ! even 
     4022               SELECT CASE(TRIM(cl_point)) 
     4023                  CASE('F','U') 
     4024                     dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0) -   & 
     4025                        &       dl_lon0(id_imin0  ,id_jmin0) ) / & 
     4026                        &     ( 2.*id_rho(jp_I) ) 
     4027                  CASE DEFAULT 
     4028                     dl_dlon=0 
     4029               END SELECT 
     4030            ELSE 
     4031               ! odd 
     4032               dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0) -   & 
     4033                  &       dl_lon0(id_imin0  ,id_jmin0) ) / & 
     4034                  &     ( 2.*id_rho(jp_I) ) 
     4035            ENDIF 
     4036 
     4037            dl_lon0F= dl_lon0(id_imin0+1,id_jmin0) + dl_dlon 
     4038            dl_lat0F= dd_lat0(id_imin0+1,id_jmin0) 
     4039 
     4040            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 
     4041            &                           dl_lon0F, dl_lat0F, 'le' ) 
     4042       
     4043            ii=il_ind(1) 
     4044 
     4045            !!!!! i-direction !!!!! 
     4046            IF( ll_even(jp_I) )THEN 
     4047               ! even 
     4048               SELECT CASE(TRIM(cl_point)) 
     4049                  CASE('T','V') 
     4050                     grid__get_fine_offset_cc(jp_I,1)=id_rho(jp_I)-ii 
     4051                  CASE DEFAULT !'F','U'  
     4052                     grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 
     4053               END SELECT 
     4054            ELSE 
     4055               ! odd 
     4056               grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 
     4057            ENDIF 
     4058 
    38564059         ELSE 
    38574060            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//& 
    3858             &                 " not match fine grid lower left corner.") 
    3859          ENDIF 
    3860          ! look for i-direction right offset 
    3861          IF( dl_lon1(il_shape1(jp_I),1) > dl_lon0(id_imax0-1,id_jmin0) )THEN 
    3862             DO ji=1,id_rho(jp_I)+2 
    3863                ii=il_shape1(jp_I)-ji+1 
    3864                IF( dl_lon1(ii,1) < dl_lon0(id_imax0-1,id_jmin0) + dp_delta )THEN 
    3865                   grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ji 
    3866                   EXIT 
    3867                ENDIF 
    3868             ENDDO 
     4061            &                 " not match fine grid left corner.") 
     4062         ENDIF 
     4063 
     4064         IF( ll_greenwich )THEN 
     4065            ! close to greenwich meridien 
     4066            ll_greenwich=.FALSE. 
     4067            ! -180:180 => 0:360 
     4068            WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0) < 0. ) 
     4069               dl_lon0(id_imin0:id_imin0+1,id_jmin0) = & 
     4070                  & dl_lon0(id_imin0:id_imin0+1,id_jmin0)+360. 
     4071            END WHERE 
     4072 
     4073            WHERE( dl_lon1(i1:i2,j1:j2) < 0. ) 
     4074               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360. 
     4075            END WHERE 
     4076         ENDIF 
     4077 
     4078         !!!!!! look for i-direction right offset !!!!!! 
     4079         i1=MAX(1,il_shape1(jp_I)-(id_rho(jp_I)+2)+1) ; i2=il_shape1(jp_I) 
     4080         j1=1                                         ; j2=1 
     4081 
     4082         ! check if cross greenwich meridien 
     4083         IF( minval(dl_lon0(id_imax0-1:id_imax0,id_jmin0))<5. .OR. & 
     4084           & maxval(dl_lon0(id_imax0-1:id_imax0,id_jmin0))>355. )THEN 
     4085            ! close to greenwich meridien 
     4086            ll_greenwich=.TRUE. 
     4087            ! 0:360 => -180:180 
     4088            WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmin0) > 180. ) 
     4089               dl_lon0(id_imax0-1:id_imax0,id_jmin0) = & 
     4090                  & dl_lon0(id_imax0-1:id_imax0,id_jmin0)-360. 
     4091            END WHERE 
     4092 
     4093            WHERE( dl_lon1(i1:i2,j1:j2) > 180. ) 
     4094               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360. 
     4095            END WHERE 
     4096         ENDIF 
     4097 
     4098         ! min lognitude of the right cell 
     4099         dl_lonmin0=dl_lon0(id_imax0-1,id_jmin0) 
     4100         IF( dl_lon1(il_shape1(jp_I),il_shape1(jp_J)) > dl_lonmin0 )THEN 
     4101 
     4102            !!!!! i-direction !!!!! 
     4103            IF( ll_even(jp_I) )THEN 
     4104               ! even 
     4105               SELECT CASE(TRIM(cl_point)) 
     4106                  CASE('F','U') 
     4107                     dl_dlon= ( dl_lon0(id_imax0  ,id_jmin0) -   & 
     4108                        &       dl_lon0(id_imax0-1,id_jmin0) ) / & 
     4109                        &     ( 2.*id_rho(jp_I) ) 
     4110                  CASE DEFAULT 
     4111                     dl_dlon=0 
     4112               END SELECT 
     4113            ELSE 
     4114               ! odd 
     4115               dl_dlon= ( dl_lon0(id_imax0  ,id_jmin0) -   & 
     4116                  &       dl_lon0(id_imax0-1,id_jmin0) ) / & 
     4117                  &     ( 2.*id_rho(jp_I) ) 
     4118            ENDIF 
     4119 
     4120            dl_lon0F= dl_lon0(id_imax0-1,id_jmin0) - dl_dlon 
     4121            dl_lat0F= dd_lat0(id_imax0-1,id_jmin0) 
     4122 
     4123            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 
     4124            &                           dl_lon0F, dl_lat0F, 'ri' ) 
     4125 
     4126            ii=(MIN(il_shape1(jp_I),(id_rho(jp_I)+2))-il_ind(1)+1) 
     4127 
     4128            !!!!! i-direction !!!!! 
     4129            IF( ll_even(jp_I) )THEN 
     4130               ! even 
     4131               SELECT CASE(TRIM(cl_point)) 
     4132                  CASE('T','V') 
     4133                     grid__get_fine_offset_cc(jp_I,2)=id_rho(jp_I)-ii 
     4134                  CASE DEFAULT !'F','U'  
     4135                     grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii 
     4136               END SELECT 
     4137            ELSE 
     4138               ! odd 
     4139               grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii 
     4140            ENDIF 
     4141 
    38694142         ELSE 
    38704143            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//& 
    3871             &                 " not match fine grid lower right corner.") 
     4144            &                 " not match fine grid right corner.") 
     4145         ENDIF 
     4146 
     4147         IF( ll_greenwich )THEN 
     4148            ! close to greenwich meridien 
     4149            ll_greenwich=.FALSE. 
     4150            ! -180:180 => 0:360 
     4151            WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmin0) < 0. ) 
     4152               dl_lon0(id_imax0-1:id_imax0,id_jmin0) = & 
     4153                  & dl_lon0(id_imax0-1:id_imax0,id_jmin0)+360. 
     4154            END WHERE 
     4155 
     4156            WHERE( dl_lon1(i1:i2,j1:j2) < 0. ) 
     4157               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360. 
     4158            END WHERE 
    38724159         ENDIF 
    38734160 
     
    38764163         grid__get_fine_offset_cc(jp_I,:)=((id_rho(jp_I)-1)/2) 
    38774164          
    3878          ! work on j-direction 
    3879  
    3880          ! look for j-direction lower offset  
    3881          IF( dd_lat1(1,1) < dd_lat0(id_imin0,id_jmin0+1) )THEN 
    3882             DO jj=1,id_rho(jp_J)+2 
    3883                IF( dd_lat1(1,jj) > dd_lat0(id_imin0,id_jmin0+1) - dp_delta )THEN 
    3884                   grid__get_fine_offset_cc(jp_J,1)=(id_rho(jp_J)+1)-jj 
    3885                   EXIT 
    3886                ENDIF 
    3887             ENDDO 
     4165         !!! work on j-direction 
     4166         !!! look for j-direction lower offset  
     4167         i1=1 ; i2=1 
     4168         j1=1 ; j2=MIN((id_rho(jp_J)+2),il_shape1(jp_J)) 
     4169 
     4170 
     4171         ! max latitude of the lower cell 
     4172         dl_latmax0=dd_lat0(id_imin0,id_jmin0+1) 
     4173         IF( dd_lat1(1,1) < dl_latmax0 )THEN 
     4174 
     4175            IF( ll_even(jp_J) )THEN 
     4176               ! even 
     4177               SELECT CASE(TRIM(cl_point)) 
     4178                  CASE('F','V') 
     4179                     dl_dlat= ( dd_lat0(id_imin0,id_jmin0+1) -   & 
     4180                        &       dd_lat0(id_imin0,id_jmin0  ) ) / & 
     4181                        &     ( 2.*id_rho(jp_J) ) 
     4182                  CASE DEFAULT 
     4183                     dl_dlat=0 
     4184               END SELECT 
     4185            ELSE 
     4186               ! odd 
     4187               dl_dlat= ( dd_lat0(id_imin0,id_jmin0+1) -   & 
     4188                  &       dd_lat0(id_imin0,id_jmin0  ) ) / & 
     4189                  &     ( 2.*id_rho(jp_J) ) 
     4190            ENDIF 
     4191 
     4192            dl_lon0F= dl_lon0(id_imin0,id_jmin0+1) 
     4193            dl_lat0F= dd_lat0(id_imin0,id_jmin0+1) + dl_dlat  
     4194             
     4195            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 
     4196            &                           dl_lon0F, dl_lat0F, 'lo' ) 
     4197 
     4198            ij=il_ind(2) 
     4199 
     4200            !!!!! i-direction !!!!! 
     4201            IF( ll_even(jp_I) )THEN 
     4202               ! even 
     4203               SELECT CASE(TRIM(cl_point)) 
     4204                  CASE('T','V') 
     4205                     grid__get_fine_offset_cc(jp_I,1)=id_rho(jp_I)-ii 
     4206                  CASE DEFAULT !'F','U'  
     4207                     grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 
     4208               END SELECT 
     4209            ELSE 
     4210               ! odd 
     4211               grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 
     4212            ENDIF 
     4213 
    38884214         ELSE 
    38894215            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//& 
    3890             &                 " not match fine grid upper left corner.") 
    3891          ENDIF 
    3892  
    3893          ! look for j-direction upper offset  
    3894          IF( dd_lat1(1,il_shape1(jp_J)) > dd_lat0(id_imin0,id_jmax0-1) )THEN 
    3895             DO jj=1,id_rho(jp_J)+2 
    3896                ij=il_shape1(jp_J)-jj+1 
    3897                IF( dd_lat1(1,ij) < dd_lat0(id_imin0,id_jmax0-1) + dp_delta )THEN 
    3898                   grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-jj 
    3899                   EXIT 
    3900                ENDIF 
    3901             ENDDO 
     4216            &                 " not match fine grid lower corner.") 
     4217         ENDIF 
     4218 
     4219         !!! look for j-direction upper offset  
     4220         i1=1                                         ; i2=1 
     4221         j1=MAX(1,il_shape1(jp_J)-(id_rho(jp_J)+2)+1) ; j2=il_shape1(jp_J) 
     4222 
     4223         ! min latitude of the upper cell 
     4224         dl_latmin0=dd_lat0(id_imin0,id_jmax0-1) 
     4225         IF( dd_lat1(il_shape1(jp_I),il_shape1(jp_J)) > dl_latmin0 )THEN 
     4226 
     4227            IF( ll_even(jp_J) )THEN 
     4228               ! even 
     4229               SELECT CASE(TRIM(cl_point)) 
     4230                  CASE('F','V') 
     4231                     dl_dlat= ( dd_lat0(id_imin0,id_jmax0  ) -   & 
     4232                        &       dd_lat0(id_imin0,id_jmax0-1) ) / & 
     4233                        &     ( 2.*id_rho(jp_J) ) 
     4234                  CASE DEFAULT 
     4235                     dl_dlat=0 
     4236               END SELECT 
     4237            ELSE 
     4238               ! odd 
     4239               dl_dlat= ( dd_lat0(id_imin0,id_jmax0  ) -   & 
     4240                  &       dd_lat0(id_imin0,id_jmax0-1) ) / & 
     4241                  &     ( 2*id_rho(jp_J) ) 
     4242            ENDIF 
     4243 
     4244            dl_lon0F= dl_lon0(id_imin0,id_jmax0-1)  
     4245            dl_lat0F= dd_lat0(id_imin0,id_jmax0-1) - dl_dlat 
     4246             
     4247            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 
     4248            &                           dl_lon0F, dl_lat0F, 'up' ) 
     4249 
     4250            ij=(MIN(il_shape1(jp_J),(id_rho(jp_J)+2))-il_ind(2)+1) 
     4251 
     4252            !!!!! j-direction !!!!! 
     4253            IF( ll_even(jp_J) )THEN 
     4254               ! even 
     4255               SELECT CASE(TRIM(cl_point)) 
     4256                  CASE('T','U') 
     4257                     grid__get_fine_offset_cc(jp_J,2)=id_rho(jp_J)-ij 
     4258                  CASE DEFAULT !'F','V' 
     4259                     grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij 
     4260               END SELECT 
     4261            ELSE 
     4262               ! odd 
     4263               grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij 
     4264            ENDIF 
     4265 
    39024266         ELSE 
    39034267            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//& 
     4268            &                 " not match fine grid upper corner.") 
     4269         ENDIF 
     4270 
     4271      ELSE ! il_shape1(1) > 1 .AND. il_shape1(2) > 1  
     4272 
     4273         !!!!!! look for lower left offset !!!!!! 
     4274         i1=1 ; i2=MIN((id_rho(jp_I)+2),il_shape1(jp_I)) 
     4275         j1=1 ; j2=MIN((id_rho(jp_J)+2),il_shape1(jp_J)) 
     4276 
     4277         ! check if cross greenwich meridien 
     4278         IF( minval(dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1))<5. .OR. & 
     4279           & maxval(dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1))>355. )THEN 
     4280            ! close to greenwich meridien 
     4281            ll_greenwich=.TRUE. 
     4282            ! 0:360 => -180:180 
     4283            WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) > 180. ) 
     4284               dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) = & 
     4285                  & dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1)-360. 
     4286            END WHERE 
     4287 
     4288            WHERE( dl_lon1(i1:i2,j1:j2) > 180. ) 
     4289               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360. 
     4290            END WHERE 
     4291         ENDIF 
     4292 
     4293         ! max longitude of the lower left cell 
     4294         dl_lonmax0=MAX(dl_lon0(id_imin0+1,id_jmin0),dl_lon0(id_imin0+1,id_jmin0+1)) 
     4295         ! max latitude of the lower left cell 
     4296         dl_latmax0=MAX(dd_lat0(id_imin0,id_jmin0+1),dd_lat0(id_imin0+1,id_jmin0+1)) 
     4297         IF( dl_lon1(1,1) < dl_lonmax0 .AND. & 
     4298           & dd_lat1(1,1) < dl_latmax0 )THEN 
     4299 
     4300            !!!!! i-direction !!!!! 
     4301            IF( ll_even(jp_I) )THEN 
     4302               ! even 
     4303               SELECT CASE(TRIM(cl_point)) 
     4304                  CASE('F','U') 
     4305                     dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0+1) -   & 
     4306                        &       dl_lon0(id_imin0  ,id_jmin0+1) ) / & 
     4307                        &     ( 2.*id_rho(jp_I) ) 
     4308                  CASE DEFAULT 
     4309                     dl_dlon=0 
     4310               END SELECT 
     4311            ELSE 
     4312               ! odd 
     4313               dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0+1) -   & 
     4314                  &       dl_lon0(id_imin0  ,id_jmin0+1) ) / & 
     4315                  &     ( 2.*id_rho(jp_I) ) 
     4316            ENDIF 
     4317 
     4318            !!!!! j-direction !!!!! 
     4319            IF( ll_even(jp_J) )THEN 
     4320               ! even 
     4321               SELECT CASE(TRIM(cl_point)) 
     4322                  CASE('F','V') 
     4323                     dl_dlat= ( dd_lat0(id_imin0+1,id_jmin0+1) -   & 
     4324                        &       dd_lat0(id_imin0+1,id_jmin0  ) ) / & 
     4325                        &     ( 2.*id_rho(jp_J) ) 
     4326                  CASE DEFAULT 
     4327                     dl_dlat=0 
     4328               END SELECT 
     4329            ELSE 
     4330               ! odd 
     4331               dl_dlat= ( dd_lat0(id_imin0+1,id_jmin0+1) -   & 
     4332                  &       dd_lat0(id_imin0+1,id_jmin0  ) ) / & 
     4333                  &     ( 2.*id_rho(jp_J) ) 
     4334            ENDIF 
     4335 
     4336            dl_lon0F= dl_lon0(id_imin0+1,id_jmin0+1) + dl_dlon 
     4337            dl_lat0F= dd_lat0(id_imin0+1,id_jmin0+1) + dl_dlat 
     4338 
     4339            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 
     4340            &                           dl_lon0F, dl_lat0F, 'll' ) 
     4341 
     4342            ii=il_ind(1) 
     4343            ij=il_ind(2) 
     4344 
     4345            !!!!! i-direction !!!!! 
     4346            IF( ll_even(jp_I) )THEN 
     4347               ! even 
     4348               SELECT CASE(TRIM(cl_point)) 
     4349                  CASE('T','V') 
     4350                     grid__get_fine_offset_cc(jp_I,1)=id_rho(jp_I)-ii 
     4351                  CASE DEFAULT !'F','U'  
     4352                     grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 
     4353               END SELECT 
     4354            ELSE 
     4355               ! odd 
     4356               grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 
     4357            ENDIF 
     4358 
     4359            !!!!! j-direction !!!!! 
     4360            IF( ll_even(jp_J) )THEN 
     4361               ! even 
     4362               SELECT CASE(TRIM(cl_point)) 
     4363                  CASE('T','U') 
     4364                     grid__get_fine_offset_cc(jp_J,1)=id_rho(jp_J)-ij 
     4365                  CASE DEFAULT !'F','V' 
     4366                     grid__get_fine_offset_cc(jp_J,1)=(id_rho(jp_J)+1)-ij 
     4367               END SELECT 
     4368            ELSE 
     4369               ! odd 
     4370               grid__get_fine_offset_cc(jp_J,1)=(id_rho(jp_J)+1)-ij 
     4371            ENDIF 
     4372 
     4373         ELSE 
     4374            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do"//& 
     4375            &                 " not match fine grid lower left corner.") 
     4376         ENDIF 
     4377 
     4378         IF( ll_greenwich )THEN 
     4379            ! close to greenwich meridien 
     4380            ll_greenwich=.FALSE. 
     4381            ! -180:180 => 0:360  
     4382            WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) < 0. ) 
     4383               dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) = & 
     4384                  & dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1)+360. 
     4385            END WHERE 
     4386 
     4387            WHERE( dl_lon1(i1:i2,j1:j2) < 0. ) 
     4388               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360. 
     4389            END WHERE 
     4390         ENDIF 
     4391 
     4392         !!!!!! look for upper right offset !!!!!! 
     4393         i1=MAX(1,il_shape1(jp_I)-(id_rho(jp_I)+2)+1) ; i2=il_shape1(jp_I) 
     4394         j1=MAX(1,il_shape1(jp_J)-(id_rho(jp_J)+2)+1) ; j2=il_shape1(jp_J) 
     4395 
     4396         ! check if cross greenwich meridien 
     4397         IF( minval(dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0))<5. .OR. & 
     4398           & maxval(dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0))>355. )THEN 
     4399            ! close to greenwich meridien 
     4400            ll_greenwich=.TRUE. 
     4401            ! 0:360 => -180:180 
     4402            WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) > 180. ) 
     4403               dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) = & 
     4404                  & dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0)-360. 
     4405            END WHERE 
     4406 
     4407            WHERE( dl_lon1(i1:i2,j1:j2) > 180. ) 
     4408               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360. 
     4409            END WHERE 
     4410         ENDIF 
     4411 
     4412         ! min latitude of the upper right cell 
     4413         dl_lonmin0=MIN(dl_lon0(id_imax0-1,id_jmax0-1),dl_lon0(id_imax0-1,id_jmax0)) 
     4414         ! min latitude of the upper right cell 
     4415         dl_latmin0=MIN(dd_lat0(id_imax0-1,id_jmax0-1),dd_lat0(id_imax0,id_jmax0-1)) 
     4416         IF( dl_lon1(il_shape1(jp_I),il_shape1(jp_J)) > dl_lonmin0 .AND. & 
     4417           & dd_lat1(il_shape1(jp_I),il_shape1(jp_J)) > dl_latmin0 )THEN 
     4418 
     4419            !!!!! i-direction !!!!! 
     4420            IF( ll_even(jp_I) )THEN 
     4421               ! even 
     4422               SELECT CASE(TRIM(cl_point)) 
     4423                  CASE('F','U') 
     4424                     dl_dlon= ( dl_lon0(id_imax0  ,id_jmax0-1) -   & 
     4425                        &       dl_lon0(id_imax0-1,id_jmax0-1) ) / & 
     4426                        &     ( 2.*id_rho(jp_I) ) 
     4427                  CASE DEFAULT 
     4428                     dl_dlon=0 
     4429               END SELECT                
     4430            ELSE 
     4431               ! odd 
     4432               dl_dlon= ( dl_lon0(id_imax0  ,id_jmax0-1) -   & 
     4433                  &       dl_lon0(id_imax0-1,id_jmax0-1) ) / & 
     4434                  &     ( 2*id_rho(jp_I) ) 
     4435            ENDIF 
     4436 
     4437            !!!!! j-direction !!!!! 
     4438            IF( ll_even(jp_J) )THEN 
     4439               ! even 
     4440               SELECT CASE(TRIM(cl_point)) 
     4441                  CASE('F','V') 
     4442                     dl_dlat= ( dd_lat0(id_imax0-1,id_jmax0  ) -   & 
     4443                        &       dd_lat0(id_imax0-1,id_jmax0-1) ) / & 
     4444                        &     ( 2.*id_rho(jp_J) ) 
     4445                  CASE DEFAULT 
     4446                     dl_dlat=0 
     4447               END SELECT 
     4448            ELSE 
     4449               ! odd 
     4450               dl_dlat= ( dd_lat0(id_imax0-1,id_jmax0  ) -   & 
     4451                  &       dd_lat0(id_imax0-1,id_jmax0-1) ) / & 
     4452                  &     ( 2*id_rho(jp_J) ) 
     4453            ENDIF 
     4454 
     4455            dl_lon0F= dl_lon0(id_imax0-1,id_jmax0-1) - dl_dlon 
     4456            dl_lat0F= dd_lat0(id_imax0-1,id_jmax0-1) - dl_dlat 
     4457 
     4458            il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 
     4459            &                           dl_lon0F, dl_lat0F, 'ur' ) 
     4460 
     4461            ii=(MIN(il_shape1(jp_I),(id_rho(jp_I)+2))-il_ind(1)+1) 
     4462            ij=(MIN(il_shape1(jp_J),(id_rho(jp_J)+2))-il_ind(2)+1) 
     4463 
     4464            !!!!! i-direction !!!!! 
     4465            IF( ll_even(jp_I) )THEN 
     4466               ! even 
     4467               SELECT CASE(TRIM(cl_point)) 
     4468                  CASE('T','V') 
     4469                     grid__get_fine_offset_cc(jp_I,2)=id_rho(jp_I)-ii 
     4470                  CASE DEFAULT !'F','U'  
     4471                     grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii 
     4472               END SELECT 
     4473            ELSE 
     4474               ! odd 
     4475               grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii 
     4476            ENDIF 
     4477 
     4478            !!!!! j-direction !!!!! 
     4479            IF( ll_even(jp_J) )THEN 
     4480               ! even 
     4481               SELECT CASE(TRIM(cl_point)) 
     4482                  CASE('T','U') 
     4483                     grid__get_fine_offset_cc(jp_J,2)=id_rho(jp_J)-ij 
     4484                  CASE DEFAULT !'F','V' 
     4485                     grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij 
     4486               END SELECT 
     4487            ELSE 
     4488               ! odd 
     4489               grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij 
     4490            ENDIF 
     4491 
     4492         ELSE 
     4493            CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do"//& 
    39044494            &                 " not match fine grid upper right corner.") 
    3905          ENDIF          
    3906  
    3907       ELSE ! il_shape1(1) > 1 .AND. il_shape1(2) > 1  
    3908  
    3909          ! look for lower left offset 
    3910          IF( dl_lon1(1,1) < dl_lon0(id_imin0+1,id_jmin0+1) )THEN 
    3911  
    3912             ii=1 
    3913             ij=1 
    3914             DO ji=1,(id_rho(jp_I)+2)*(id_rho(jp_J)+2) 
    3915  
    3916                ll_ii=.FALSE. 
    3917                ll_ij=.FALSE. 
    3918  
    3919                IF( dl_lon1(ii,ij) >= dl_lon0(id_imin0+1,id_jmin0+1)-dp_delta .AND. & 
    3920                &   dd_lat1(ii,ij) >= dd_lat0(id_imin0+1,id_jmin0+1)-dp_delta )THEN 
    3921                   grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 
    3922                   grid__get_fine_offset_cc(jp_J,1)=(id_rho(jp_J)+1)-ij 
    3923                   EXIT 
    3924                ENDIF 
    3925  
    3926                IF( dl_lon1(ii+1,ij) <= dl_lon0(id_imin0+1,id_jmin0+1)+dp_delta .AND. & 
    3927                &   dd_lat1(ii+1,ij) <= dd_lat0(id_imin0+1,id_jmin0+1)+dp_delta )THEN 
    3928                   ll_ii=.TRUE. 
    3929                ENDIF 
    3930                IF( dl_lon1(ii,ij+1) <= dl_lon0(id_imin0+1,id_jmin0+1)+dp_delta .AND. & 
    3931                &   dd_lat1(ii,ij+1) <= dd_lat0(id_imin0+1,id_jmin0+1)+dp_delta )THEN 
    3932                   ll_ij=.TRUE. 
    3933                ENDIF 
    3934  
    3935                IF( ll_ii ) ii=ii+1 
    3936                IF( ll_ij ) ij=ij+1 
    3937  
    3938             ENDDO 
    3939  
    3940          ELSE 
    3941             CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//& 
    3942             &                 " not match fine grid lower left corner.") 
    3943          ENDIF 
    3944  
    3945          ! look for upper right offset 
    3946          IF( dl_lon1(il_shape1(jp_I),il_shape1(jp_J)) > & 
    3947             & dl_lon0(id_imax0-1,id_jmax0-1) )THEN 
    3948  
    3949             ii=il_shape1(jp_I) 
    3950             ij=il_shape1(jp_J) 
    3951             DO ji=1,(id_rho(jp_I)+2)*(id_rho(jp_J)+2) 
    3952  
    3953                ll_ii=.FALSE. 
    3954                ll_ij=.FALSE. 
    3955  
    3956                IF( dl_lon1(ii,ij) <= dl_lon0(id_imax0-1,id_jmax0-1)+dp_delta .AND. & 
    3957                &   dd_lat1(ii,ij) <= dd_lat0(id_imax0-1,id_jmax0-1)+dp_delta )THEN 
    3958                   grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-(il_shape1(jp_I)+1-ii) 
    3959                   grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-(il_shape1(jp_J)+1-ij) 
    3960                   EXIT 
    3961                ENDIF 
    3962  
    3963                IF( dl_lon1(ii-1,ij) >= dl_lon0(id_imax0-1,id_jmax0-1)-dp_delta .AND. & 
    3964                &   dd_lat1(ii-1,ij) >= dd_lat0(id_imax0-1,id_jmax0-1)-dp_delta )THEN 
    3965                   ll_ii=.TRUE. 
    3966                ENDIF 
    3967                IF( dl_lon1(ii,ij-1) >= dl_lon0(id_imax0-1,id_jmax0-1)-dp_delta .AND. & 
    3968                &   dd_lat1(ii,ij-1) >= dd_lat0(id_imax0-1,id_jmax0-1)-dp_delta )THEN 
    3969                   ll_ij=.TRUE. 
    3970                ENDIF 
    3971  
    3972                IF( ll_ii ) ii=ii-1 
    3973                IF( ll_ij ) ij=ij-1 
    3974  
    3975             ENDDO 
    3976  
    3977          ELSE 
    3978             CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//& 
    3979             &                 " not match fine grid upper right corner.") 
     4495         ENDIF 
     4496 
     4497         IF( ll_greenwich )THEN 
     4498            ! close to greenwich meridien 
     4499            ll_greenwich=.FALSE. 
     4500            ! -180:180 => 0:360 
     4501            WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) < 0. ) 
     4502               dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) = & 
     4503                  & dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0)+360. 
     4504            END WHERE 
     4505 
     4506            WHERE( dl_lon1(i1:i2,j1:j2) < 0. ) 
     4507               dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360. 
     4508            END WHERE 
    39804509         ENDIF 
    39814510 
     
    39844513      DEALLOCATE( dl_lon0 ) 
    39854514      DEALLOCATE( dl_lon1 ) 
     4515 
     4516      IF( ANY(grid__get_fine_offset_cc(:,:)==-1) )THEN 
     4517         CALL logger_fatal("GRID GET FINE OFFSET: can not found "//& 
     4518         &                 " offset between coarse and fine grid.") 
     4519      ENDIF 
    39864520 
    39874521   END FUNCTION grid__get_fine_offset_cc 
     
    39954529   !> @date October, 2014 
    39964530   !> - work on mpp file structure instead of file structure 
    3997    ! 
     4531   !> @date February, 2016 
     4532   !> - use F-point to check coincidence for even refinment 
     4533   !> - use F-point estimation, if can not read it. 
     4534   !> 
    39984535   !> @param[in] td_coord0 coarse grid coordinate file structure  
    39994536   !> @param[in] td_coord1 fine   grid coordinate file structure  
     
    40204557 
    40214558      ! local variable 
    4022       INTEGER(i4) :: il_imid1 
    4023       INTEGER(i4) :: il_jmid1 
     4559      INTEGER(i4)               :: il_imid1 
     4560      INTEGER(i4)               :: il_jmid1 
    40244561       
    4025       INTEGER(i4) :: il_ew0 
    4026       INTEGER(i4) :: il_ew1 
    4027  
    4028       INTEGER(i4) :: il_imin1 
    4029       INTEGER(i4) :: il_imax1 
    4030       INTEGER(i4) :: il_jmin1 
    4031       INTEGER(i4) :: il_jmax1 
    4032  
    4033       INTEGER(i4), DIMENSION(2) :: il_indC 
    4034       INTEGER(i4), DIMENSION(2) :: il_indF 
    4035       INTEGER(i4), DIMENSION(2) :: il_iind 
    4036       INTEGER(i4), DIMENSION(2) :: il_jind 
    4037  
    4038       REAL(dp)    :: dl_lon0 
    4039       REAL(dp)    :: dl_lat0 
    4040       REAL(dp)    :: dl_lon1 
    4041       REAL(dp)    :: dl_lat1 
    4042  
    4043       REAL(dp)    :: dl_lon1p 
    4044       REAL(dp)    :: dl_lat1p 
    4045  
    4046       LOGICAL     :: ll_coincidence 
    4047  
    4048       TYPE(TVAR)  :: tl_lon0 
    4049       TYPE(TVAR)  :: tl_lat0 
    4050       TYPE(TVAR)  :: tl_lon1 
    4051       TYPE(TVAR)  :: tl_lat1 
    4052  
    4053       TYPE(TMPP)  :: tl_coord0 
    4054       TYPE(TMPP)  :: tl_coord1 
    4055  
    4056       TYPE(TDOM)  :: tl_dom0 
     4562      INTEGER(i4)               :: il_ew0 
     4563      INTEGER(i4)               :: il_ew1 
     4564 
     4565      INTEGER(i4)               :: il_ind 
     4566 
     4567      INTEGER(i4)               :: il_imin1 
     4568      INTEGER(i4)               :: il_imax1 
     4569      INTEGER(i4)               :: il_jmin1 
     4570      INTEGER(i4)               :: il_jmax1 
     4571 
     4572      INTEGER(i4), DIMENSION(2) :: il_ind0 
     4573      INTEGER(i4), DIMENSION(2) :: il_ind1 
     4574 
     4575      INTEGER(i4), DIMENSION(2) :: il_ill1 
     4576      INTEGER(i4), DIMENSION(2) :: il_ilr1 
     4577      INTEGER(i4), DIMENSION(2) :: il_iul1 
     4578      INTEGER(i4), DIMENSION(2) :: il_iur1 
     4579 
     4580      REAL(dp)                  :: dl_lon0F 
     4581      REAL(dp)                  :: dl_lat0F 
     4582      REAL(dp)                  :: dl_lon0 
     4583      REAL(dp)                  :: dl_lat0 
     4584      REAL(dp)                  :: dl_lon1F 
     4585      REAL(dp)                  :: dl_lat1F 
     4586      REAL(dp)                  :: dl_lon1 
     4587      REAL(dp)                  :: dl_lat1 
     4588 
     4589      REAL(dp)                  :: dl_delta 
     4590 
     4591      LOGICAL                   :: ll_coincidence 
     4592      LOGICAL                   :: ll_even 
     4593      LOGICAL                   :: ll_grid0F 
     4594      LOGICAL                   :: ll_grid1F 
     4595 
     4596      TYPE(TVAR)                :: tl_lon0 
     4597      TYPE(TVAR)                :: tl_lat0 
     4598      TYPE(TVAR)                :: tl_lon0F 
     4599      TYPE(TVAR)                :: tl_lat0F 
     4600      TYPE(TVAR)                :: tl_lon1 
     4601      TYPE(TVAR)                :: tl_lat1 
     4602      TYPE(TVAR)                :: tl_lon1F 
     4603      TYPE(TVAR)                :: tl_lat1F 
     4604 
     4605      TYPE(TMPP)                :: tl_coord0 
     4606      TYPE(TMPP)                :: tl_coord1 
     4607 
     4608      TYPE(TDOM)                :: tl_dom0 
    40574609 
    40584610      ! loop indices 
     
    40634615      ll_coincidence=.TRUE. 
    40644616 
     4617      ll_even=.FALSE. 
     4618      IF( MOD(id_rho(jp_I)*id_rho(jp_J),2) == 0 )THEN 
     4619         ll_even=.TRUE. 
     4620      ENDIF 
     4621 
    40654622      ! copy structure 
    40664623      tl_coord0=mpp_copy(td_coord0) 
     
    40754632 
    40764633      ! read variable value on domain 
    4077       tl_lon0=iom_dom_read_var(tl_coord0,'longitude',tl_dom0) 
    4078       tl_lat0=iom_dom_read_var(tl_coord0,'latitude' ,tl_dom0) 
     4634      il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'longitude_T') 
     4635      IF( il_ind /= 0 )THEN 
     4636         tl_lon0=iom_dom_read_var(tl_coord0,'longitude_T',tl_dom0) 
     4637      ELSE 
     4638         tl_lon0=iom_dom_read_var(tl_coord0,'longitude',tl_dom0) 
     4639      ENDIF 
     4640 
     4641      il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'latitude_T') 
     4642      IF( il_ind /= 0 )THEN 
     4643         tl_lat0=iom_dom_read_var(tl_coord0,'latitude_T' ,tl_dom0) 
     4644      ELSE 
     4645         tl_lat0=iom_dom_read_var(tl_coord0,'latitude' ,tl_dom0) 
     4646      ENDIF 
     4647 
     4648      IF( ll_even )THEN 
     4649         ! look for variable value on domain for F point 
     4650         il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'longitude_F') 
     4651         IF( il_ind /= 0 )THEN 
     4652            tl_lon0F=iom_dom_read_var(tl_coord0,'longitude_F',tl_dom0) 
     4653         ENDIF 
     4654 
     4655         il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'latitude_F') 
     4656         IF( il_ind /= 0 )THEN 
     4657            tl_lat0F=iom_dom_read_var(tl_coord0,'latitude_F' ,tl_dom0) 
     4658         ENDIF 
     4659 
     4660         ll_grid0F=.FALSE. 
     4661         IF( ASSOCIATED(tl_lon0F%d_value) .AND. & 
     4662         &   ASSOCIATED(tl_lat0F%d_value) )THEN 
     4663            ll_grid0F=.TRUE. 
     4664         ENDIF 
     4665 
     4666      ENDIF 
    40794667 
    40804668      ! close mpp files 
     
    40924680 
    40934681      ! read fine longitue and latitude 
    4094       tl_lon1=iom_mpp_read_var(tl_coord1,'longitude') 
    4095       tl_lat1=iom_mpp_read_var(tl_coord1,'latitude') 
     4682      il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), TRIM(tl_lon0%c_longname)) 
     4683      IF( il_ind /= 0 )THEN 
     4684         tl_lon1=iom_mpp_read_var(tl_coord1,TRIM(tl_lon0%c_longname)) 
     4685      ELSE 
     4686         tl_lon1=iom_mpp_read_var(tl_coord1,'longitude') 
     4687      ENDIF 
     4688      il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), TRIM(tl_lat0%c_longname)) 
     4689      IF( il_ind /= 0 )THEN 
     4690         tl_lat1=iom_mpp_read_var(tl_coord1,TRIM(tl_lat0%c_longname)) 
     4691      ELSE 
     4692         tl_lat1=iom_mpp_read_var(tl_coord1,'latitude') 
     4693      ENDIF 
    40964694       
     4695      IF( ll_even )THEN 
     4696 
     4697         ! look for variable value on domain for F point 
     4698         il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), 'longitude_F') 
     4699         IF( il_ind /= 0 )THEN 
     4700            tl_lon1F=iom_mpp_read_var(tl_coord1,'longitude_F') 
     4701         ENDIF 
     4702 
     4703         il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), 'latitude_F') 
     4704         IF( il_ind /= 0 )THEN 
     4705            tl_lat1F=iom_mpp_read_var(tl_coord1,'latitude_F') 
     4706         ENDIF 
     4707 
     4708         ll_grid1F=.FALSE. 
     4709         IF( ASSOCIATED(tl_lon1F%d_value) .AND. & 
     4710         &   ASSOCIATED(tl_lat1F%d_value) )THEN 
     4711            ll_grid1F=.TRUE. 
     4712         ENDIF 
     4713 
     4714      ENDIF 
     4715 
    40974716      ! close mpp files 
    4098       CALL iom_dom_close(tl_coord1) 
     4717      CALL iom_mpp_close(tl_coord1) 
    40994718      ! clean structure 
    41004719      CALL mpp_clean(tl_coord1) 
     
    41584777         IF( .NOT. ll_coincidence )THEN 
    41594778            CALL logger_fatal("GRID CHECK COINCIDENCE: no coincidence "//& 
    4160             &              "between fine grid and coarse grid. invalid domain" ) 
     4779            &              "between fine grid and coarse grid: invalid domain." ) 
    41614780         ENDIF 
    41624781 
     
    41724791 
    41734792      ! select closest point on coarse grid 
    4174       il_indC(:)=grid_get_closest(tl_lon0%d_value(:,:,1,1),& 
     4793      il_ind0(:)=grid_get_closest(tl_lon0%d_value(:,:,1,1),& 
    41754794      &                           tl_lat0%d_value(:,:,1,1),& 
    41764795      &                           dl_lon1, dl_lat1   ) 
    41774796 
    4178       IF( ANY(il_indC(:)==0) )THEN 
     4797      IF( ANY(il_ind0(:)==0) )THEN 
    41794798         CALL logger_fatal("GRID CHECK COINCIDENCE: can not find valid "//& 
    4180          &              "coarse grid indices. invalid domain" ) 
    4181       ENDIF 
    4182  
    4183       dl_lon0=tl_lon0%d_value(il_indC(1),il_indC(2),1,1) 
    4184       dl_lat0=tl_lat0%d_value(il_indC(1),il_indC(2),1,1) 
    4185  
    4186       ! look for closest fine grid point from selected coarse grid point 
    4187       il_iind(:)=MAXLOC( tl_lon1%d_value(:,:,1,1), & 
    4188       &                  tl_lon1%d_value(:,:,1,1) <= dl_lon0 ) 
    4189  
    4190       il_jind(:)=MAXLOC( tl_lat1%d_value(:,:,1,1), & 
    4191       &                  tl_lat1%d_value(:,:,1,1) <= dl_lat0 ) 
    4192  
    4193       il_indF(1)=il_iind(1) 
    4194       il_indF(2)=il_jind(2) 
    4195  
    4196       IF( ANY(il_indF(:)==0) )THEN 
    4197          CALL logger_fatal("GRID CHECK COINCIDENCE: can not find valid "//& 
    4198          &              "fine grid indices. invalid domain" ) 
    4199       ENDIF 
    4200  
    4201       dl_lon1=tl_lon1%d_value(il_indF(1),il_indF(2),1,1) 
    4202       dl_lat1=tl_lat1%d_value(il_indF(1),il_indF(2),1