Changeset 9023


Ignore:
Timestamp:
2017-12-13T18:08:50+01:00 (3 years ago)
Author:
timgraham
Message:

Merged METO_MERCATOR branch and resolved all conflicts in OPA_SRC

Location:
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC
Files:
47 edited
1 copied

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r9019 r9023  
    2020   !!   dyn_asm_inc    : Apply the dynamic (u and v) increments 
    2121   !!   ssh_asm_inc    : Apply the SSH increment 
     22   !!   ssh_asm_div    : Apply divergence associated with SSH increment 
    2223   !!   seaice_asm_inc : Apply the seaice increment 
    2324   !!---------------------------------------------------------------------- 
     
    4849   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments 
    4950   PUBLIC   ssh_asm_inc    !: Apply the SSH increment 
     51   PUBLIC   ssh_asm_div    !: Apply the SSH divergence 
    5052   PUBLIC   seaice_asm_inc !: Apply the seaice increment 
    5153 
     
    785787   END SUBROUTINE ssh_asm_inc 
    786788 
     789   SUBROUTINE ssh_asm_div( kt, phdivn ) 
     790      !!---------------------------------------------------------------------- 
     791      !!                  ***  ROUTINE ssh_asm_div  *** 
     792      !! 
     793      !! ** Purpose :   ssh increment with z* is incorporated via a correction of the local divergence           
     794      !!                across all the water column 
     795      !! 
     796      !! ** Method  : 
     797      !!                CAUTION : sshiau is positive (inflow) decreasing the 
     798      !!                          divergence and expressed in m/s 
     799      !! 
     800      !! ** Action  :   phdivn   decreased by the ssh increment 
     801      !!---------------------------------------------------------------------- 
     802      INTEGER, INTENT(IN) :: kt                               ! ocean time-step index 
     803      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence 
     804      !! 
     805      INTEGER  ::   jk                                        ! dummy loop index 
     806      REAL(wp), DIMENSION(:,:)  , POINTER       ::   ztim     ! local array 
     807      !!---------------------------------------------------------------------- 
     808      !  
     809#if defined key_asminc 
     810      CALL ssh_asm_inc( kt ) !==   (calculate increments) 
     811      ! 
     812      IF( ln_linssh ) THEN  
     813         phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t_n(:,:,1) * tmask(:,:,1) 
     814      ELSE  
     815         CALL wrk_alloc( jpi,jpj, ztim) 
     816         ztim(:,:) = ssh_iau(:,:) / ( ht_n(:,:) + 1.0 - ssmask(:,:) ) 
     817         DO jk = 1, jpkm1                                  
     818            phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk)  
     819         END DO 
     820         ! 
     821         CALL wrk_dealloc( jpi,jpj, ztim) 
     822      ENDIF 
     823#endif 
     824      ! 
     825   END SUBROUTINE ssh_asm_div 
    787826 
    788827   SUBROUTINE seaice_asm_inc( kt, kindic ) 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90

    r9019 r9023  
    2121   USE phycst          ! physical constants 
    2222   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     23   USE wet_dry         ! Use wet dry to get reference ssh level 
    2324   USE in_out_manager  ! 
    2425 
     
    169170         ii = idx%nbi(jb,igrd) 
    170171         ij = idx%nbj(jb,igrd) 
    171          spgu(ii, ij) = dta%ssh(jb) 
     172         IF( ll_wd ) THEN 
     173            spgu(ii, ij) = dta%ssh(jb)  - ssh_ref  
     174         ELSE 
     175            spgu(ii, ij) = dta%ssh(jb) 
     176         ENDIF 
    172177      END DO 
    173178 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r9019 r9023  
    146146         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    147147 
    148       CALL iom_put( "ssh" , sshn )                 ! sea surface height 
     148      IF( ll_wd ) THEN 
     149         CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying) 
     150      ELSE 
     151         CALL iom_put( "ssh" , sshn )              ! sea surface height 
     152      ENDIF 
     153 
    149154      IF( iom_use("wetdep") )   &                  ! wet depth 
    150          CALL iom_put( "wetdep" , ht_wd(:,:) + sshn(:,:) ) 
     155         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) ) 
    151156       
    152157      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r9019 r9023  
    5353   INTEGER,  PUBLIC :: nn_baro          !: Number of barotropic iterations during one baroclinic step (rdt) 
    5454   REAL(wp), PUBLIC :: rn_bt_cmax       !: Maximum allowed courant number (used if ln_bt_auto=T) 
     55   REAL(wp), PUBLIC :: rn_bt_alpha      !: Time stepping diffusion parameter 
    5556 
    5657 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r9019 r9023  
    3838   USE c1d            ! 1D configuration 
    3939   USE dyncor_c1d     ! 1D configuration: Coriolis term    (cor_c1d routine) 
    40    USE wet_dry        ! wetting and drying 
     40   USE wet_dry,  ONLY : ll_wd 
    4141   ! 
    4242   USE in_out_manager ! I/O manager 
     
    667667      ENDIF 
    668668      ! 
    669       IF( ln_wd ) THEN              ! wetting and drying domain 
     669      IF( ll_wd ) THEN              ! wetting and drying domain 
    670670         CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 ) 
    671          CALL iom_rstput( 0, 0, inum, 'ht_wd'  , ht_wd  , ktype = jp_r8 ) 
    672671      ENDIF 
    673672      ! 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r9019 r9023  
    275275               ENDIF 
    276276            END DO 
     277#if defined key_agrif  
     278       IF( .NOT. AGRIF_Root() ) THEN  
     279          IF ((nbondi ==  1).OR.(nbondi == 2)) fmask(nlci-1 , :     ,jk) = 0.e0      ! east  
     280          IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1      , :     ,jk) = 0.e0      ! west  
     281          IF ((nbondj ==  1).OR.(nbondj == 2)) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north  
     282          IF ((nbondj == -1).OR.(nbondj == 2)) fmask(:      ,1      ,jk) = 0.e0      ! south  
     283       ENDIF  
     284#endif  
    277285         END DO 
    278286         ! 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r9019 r9023  
    672672      ! 
    673673      INTEGER ::   ji, jj, jk                                       ! dummy loop indices 
    674       REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd = T/F 
     674      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F 
    675675      !!---------------------------------------------------------------------- 
    676676      ! 
    677677      IF( ln_timing )   CALL timing_start('dom_vvl_interpol') 
    678678      ! 
    679       IF(ln_wd) THEN 
     679      IF(ln_wd_il) THEN 
    680680        zlnwd = 1.0_wp 
    681681      ELSE 
     
    869869         ELSE                                   !* Initialize at "rest" 
    870870            ! 
    871             IF( ln_wd .AND. ( cn_cfg == 'wad' ) ) THEN 
    872               ! Wetting and drying test case 
    873               CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  ) 
    874                        tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones 
    875                        sshn (:,:)     = sshb(:,:) 
    876                        un   (:,:,:)   = ub  (:,:,:) 
    877                        vn   (:,:,:)   = vb  (:,:,:) 
    878                                                 ! uniform T-S fields and initial ssh slope 
    879                ! needs to be called here and in istate which is called later. 
    880                ! Adjust vertical metrics 
     871 
     872            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential  
     873               ! 
     874               IF( cn_cfg == 'wad' ) THEN 
     875                  ! Wetting and drying test case 
     876                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  ) 
     877                  tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones 
     878                  sshn (:,:)     = sshb(:,:) 
     879                  un   (:,:,:)   = ub  (:,:,:) 
     880                  vn   (:,:,:)   = vb  (:,:,:) 
     881               ELSE 
     882                  ! if not test case 
     883                  sshn(:,:) = -ssh_ref 
     884                  sshb(:,:) = -ssh_ref 
     885 
     886                  DO jj = 1, jpj 
     887                     DO ji = 1, jpi 
     888                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
     889 
     890                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
     891                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
     892                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
     893                        ENDIF 
     894                     ENDDO 
     895                  ENDDO 
     896               ENDIF !If test case else 
     897 
     898               ! Adjust vertical metrics for all wad 
    881899               DO jk = 1, jpk 
    882                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 
     900                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) & 
    883901                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    884                     &            + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk)) 
     902                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 
    885903               END DO 
    886904               e3t_b(:,:,:) = e3t_n(:,:,:) 
    887                ! 
    888             ELSEIF( ln_wd ) THEN 
    889                ! 
    890               DO jj = 1, jpj 
    891                 DO ji = 1, jpi 
    892                   IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN 
    893                      ! potential bug 
    894                      ! Warning this assumes 2 layers only over wetting locations. needs investigating 
    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 
    898                      sshb(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)           !!gm I don't understand that ! 
    899                      sshn(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 
    900                      ssha(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 
    901                   ENDIF 
    902                 ENDDO 
    903               ENDDO 
     905 
     906               DO ji = 1, jpi 
     907                  DO jj = 1, jpj 
     908                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
     909                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
     910                     ENDIF 
     911                  END DO  
     912               END DO  
     913  
    904914               ! 
    905915            ELSE 
     
    909919               sshn(:,:) = 0.0_wp 
    910920               ! 
    911             END IF 
     921            END IF           ! end of ll_wd edits 
    912922 
    913923            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 
     
    10221032      ! 
    10231033#if defined key_agrif 
    1024       IF(.NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' ) 
     1034      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' ) 
    10251035#endif 
    10261036      ! 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r9019 r9023  
    1818   USE dom_oce         ! ocean space and time domain 
    1919   USE phycst ,   ONLY :   rsmall 
    20    USE wet_dry,   ONLY :   ln_wd, ht_wd 
     20   USE wet_dry,   ONLY :   ll_wd  ! Wetting and drying 
    2121   ! 
    2222   USE in_out_manager  ! I/O manager 
     
    194194      ENDIF 
    195195      ! 
    196       IF( ln_wd ) THEN                                          ! wetting and drying domain 
    197          CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 ) 
    198          CALL iom_rstput( 0, 0, inum, 'ht_wd'  , ht_wd  , ktype = jp_r8 ) 
    199       ENDIF 
     196      IF( ll_wd ) CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 ) 
     197 
    200198      !                                     ! ============================ 
    201199      CALL iom_close( inum )                !        close the files  
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r9019 r9023  
    3030   USE usrdef_zgr     ! user defined vertical coordinate system 
    3131   USE depth_e3       ! depth <=> e3 
    32    USE wet_dry, ONLY: ln_wd, ht_wd 
     32   USE wet_dry,   ONLY: ll_wd, ssh_ref  ! Wetting and drying 
    3333   ! 
    3434   USE in_out_manager ! I/O manager 
     
    257257      k_bot(:,:) = INT( z2d(:,:) ) 
    258258      ! 
    259       ! bathymetry with orography (wetting and drying only) 
    260       IF( ln_wd )  CALL iom_get( inum, jpdom_data, 'ht_wd' , ht_wd  , lrowattr=ln_use_jattr ) 
     259      ! reference depth for negative bathy (wetting and drying only) 
     260      IF( ll_wd )  CALL iom_get( inum,  'rn_wd_ref_depth' , ssh_ref  ) 
    261261      ! 
    262262      CALL iom_close( inum ) 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90

    r9019 r9023  
    2525   USE iscplhsb        ! ice sheet / ocean coupling 
    2626   USE iscplini        ! ice sheet / ocean coupling 
     27#if defined key_asminc    
     28   USE asminc          ! Assimilation increment 
     29#endif 
    2730   ! 
    2831   USE in_out_manager  ! I/O manager 
     
    9396      IF( ln_rnf )   CALL sbc_rnf_div( hdivn )              !==  runoffs    ==!   (update hdivn field) 
    9497      ! 
    95       IF( ln_isf )   CALL sbc_isf_div( hdivn )              !==  ice shelf  ==!   (update hdivn field) 
     98#if defined key_asminc  
     99      IF( ln_sshinc .AND. ln_asmiau )   CALL ssh_asm_div( kt, hdivn )   !==  SSH assimilation  ==!   (update hdivn field) 
     100      !  
     101#endif 
     102      IF( ln_isf )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field) 
    96103      ! 
    97104      IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( hdivn ) !==  ice sheet  ==!   (update hdivn field) 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r9019 r9023  
    422422      !!---------------------------------------------------------------------- 
    423423      ! 
     424      IF( ln_wd_il ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
     425      ! 
    424426      IF( kt == nit000 ) THEN 
    425427         IF(lwp) WRITE(numout,*) 
     
    433435      ENDIF 
    434436      ! 
    435       IF( ln_wd ) THEN 
    436          ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) 
    437          DO jj = 2, jpjm1 
    438             DO ji = 2, jpim1  
    439                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    440                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            & 
    441                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 
    442                   &                                                         > rn_wdmin1 + rn_wdmin2 
    443                ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             & 
    444                   &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    445                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    446  
    447                IF(ll_tmp1) THEN 
    448                   zcpx(ji,jj) = 1.0_wp 
    449                ELSE IF(ll_tmp2) THEN 
    450                   ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    451                   zcpx(ji,jj) = ABS(   ( sshn(ji+1,jj)+ht_wd(ji+1,jj) - sshn(ji,jj)-ht_wd(ji,jj) )   & 
    452                      &                / ( sshn(ji+1,jj)                - sshn(ji,jj)              )   ) 
    453                ELSE 
    454                   zcpx(ji,jj) = 0._wp 
    455                ENDIF 
    456                ! 
    457                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    458                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            & 
    459                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 
    460                   &                                                         > rn_wdmin1 + rn_wdmin2 
    461                ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             & 
    462                   &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    463                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     437      IF( ln_wd_il ) THEN 
     438        DO jj = 2, jpjm1 
     439           DO ji = 2, jpim1  
     440             ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                & 
     441                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
     442                  &    MAX(  sshn(ji,jj) +  ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     443                  &                                                       > rn_wdmin1 + rn_wdmin2 
     444             ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (       & 
     445                  &    MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     446                  &    MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     447 
     448             IF(ll_tmp1) THEN 
     449               zcpx(ji,jj) = 1.0_wp 
     450             ELSE IF(ll_tmp2) THEN 
     451               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     452               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
     453                           &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
     454             ELSE 
     455               zcpx(ji,jj) = 0._wp 
     456             END IF 
     457       
     458             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     459                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
     460                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
     461                  &                                                      > rn_wdmin1 + rn_wdmin2 
     462             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (        & 
     463                  &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                & 
     464                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    464465                  ! 
    465466               IF(ll_tmp1) THEN 
     
    476477         CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    477478      ENDIF 
     479 
     480             IF(ll_tmp1) THEN 
     481               zcpy(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) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
     485                           &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     486             ELSE 
     487               zcpy(ji,jj) = 0._wp 
     488             END IF 
     489           END DO 
     490        END DO 
     491        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     492      END IF 
    478493 
    479494      ! Surface value 
     
    491506               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 
    492507            ! 
    493             IF( ln_wd ) THEN 
     508            IF( ln_wd_il ) THEN 
    494509               zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
    495510               zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     
    521536                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 
    522537               ! 
    523                IF( ln_wd ) THEN 
     538               IF( ln_wd_il ) THEN 
    524539                  zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
    525540                  zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     
    535550      END DO 
    536551      ! 
    537       IF( ln_wd )   DEALLOCATE( zcpx , zcpy ) 
     552      IF( ln_wd_il )   DEALLOCATE( zcpx , zcpy ) 
    538553      ! 
    539554   END SUBROUTINE hpg_sco 
     
    667682      !!---------------------------------------------------------------------- 
    668683      ! 
    669       IF( ln_wd ) THEN 
     684      IF( ln_wd_il ) THEN 
    670685         ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) 
    671          DO jj = 2, jpjm1 
    672             DO ji = 2, jpim1  
    673                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    674                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            & 
    675                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 
    676                   &                                                         > rn_wdmin1 + rn_wdmin2 
    677                ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             & 
    678                   &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    679                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    680  
    681                IF(ll_tmp1) THEN 
    682                   zcpx(ji,jj) = 1.0_wp 
    683                ELSE IF(ll_tmp2) THEN 
    684                   ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    685                   zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 
    686                               &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
    687                ELSE 
    688                   zcpx(ji,jj) = 0._wp 
    689                ENDIF 
     686        DO jj = 2, jpjm1 
     687           DO ji = 2, jpim1  
     688             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     689                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            & 
     690                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     691                  &                                                      > rn_wdmin1 + rn_wdmin2 
     692             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (        & 
     693                  &    MAX(   sshn(ji,jj)             ,  sshn(ji+1,jj) ) >                & 
     694                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     695             IF(ll_tmp1) THEN 
     696               zcpx(ji,jj) = 1.0_wp 
     697             ELSE IF(ll_tmp2) THEN 
     698               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     699               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
     700                           &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
     701             ELSE 
     702               zcpx(ji,jj) = 0._wp 
     703             END IF 
    690704       
    691                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    692                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            & 
    693                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 
    694                   &                                                         > rn_wdmin1 + rn_wdmin2 
    695                ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             & 
    696                   &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    697                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    698  
    699                IF(ll_tmp1) THEN 
    700                   zcpy(ji,jj) = 1.0_wp 
    701                ELSE IF(ll_tmp2) THEN 
    702                   ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    703                   zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 
    704                               &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
    705                ELSE 
    706                   zcpy(ji,jj) = 0._wp 
    707                ENDIF 
    708             END DO 
    709          END DO 
    710          CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    711       ENDIF 
     705             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     706                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
     707                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 
     708                  &                                                      > rn_wdmin1 + rn_wdmin2 
     709             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (        & 
     710                  &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                & 
     711                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     712 
     713             IF(ll_tmp1) THEN 
     714               zcpy(ji,jj) = 1.0_wp 
     715             ELSE IF(ll_tmp2) THEN 
     716               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
     717               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
     718                           &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     719             ELSE 
     720               zcpy(ji,jj) = 0._wp 
     721             END IF 
     722           END DO 
     723        END DO 
     724        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     725      END IF 
    712726 
    713727      IF( kt == nit000 ) THEN 
     
    880894            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 
    881895            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 
    882             IF( ln_wd ) THEN 
     896            IF( ln_wd_il ) THEN 
    883897              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
    884898              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     
    903917                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    904918                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj) 
    905                IF( ln_wd ) THEN 
     919               IF( ln_wd_il ) THEN 
    906920                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
    907921                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     
    914928      END DO 
    915929      ! 
    916       IF( ln_wd )   DEALLOCATE( zcpx, zcpy ) 
     930      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy ) 
    917931      ! 
    918932   END SUBROUTINE hpg_djc 
     
    959973      IF( ln_linssh )   znad = 0._wp 
    960974 
    961       IF( ln_wd ) THEN 
     975      IF( ln_wd_il ) THEN 
    962976         ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) 
    963977         DO jj = 2, jpjm1 
    964             DO ji = 2, jpim1  
    965                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    966                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            & 
    967                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 
    968                   &                                                         > rn_wdmin1 + rn_wdmin2 
    969                ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             & 
    970                   &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    971                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    972  
    973                IF(ll_tmp1) THEN 
    974                   zcpx(ji,jj) = 1.0_wp 
    975                ELSE IF(ll_tmp2) THEN 
    976                   ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    977                   zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 
    978                              &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
    979                ELSE 
    980                   zcpx(ji,jj) = 0._wp 
    981                ENDIF 
     978           DO ji = 2, jpim1  
     979             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     980                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            & 
     981                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     982                  &                                                      > rn_wdmin1 + rn_wdmin2 
     983             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (         & 
     984                  &    MAX(   sshn(ji,jj)             ,  sshn(ji+1,jj) ) >                & 
     985                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     986 
     987             IF(ll_tmp1) THEN 
     988               zcpx(ji,jj) = 1.0_wp 
     989             ELSE IF(ll_tmp2) THEN 
     990               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     991               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
     992                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
     993               
     994                zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
     995             ELSE 
     996               zcpx(ji,jj) = 0._wp 
     997             END IF 
    982998       
    983                ll_tmp1 = MIN(   sshn(ji,jj)             ,   sshn(ji,jj+1) ) >                & 
    984                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            & 
    985                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 
    986                   &                                                         > rn_wdmin1 + rn_wdmin2 
    987                ll_tmp2 = ( ABS( sshn(ji,jj)             -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             & 
    988                   &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    989                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    990  
    991                IF(ll_tmp1) THEN 
    992                   zcpy(ji,jj) = 1.0_wp 
    993                ELSE IF(ll_tmp2) THEN 
    994                   ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    995                   zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 
    996                            &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) ) 
     999             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     1000                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
     1001                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
     1002                  &                                                      > rn_wdmin1 + rn_wdmin2 
     1003             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (      & 
     1004                  &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                & 
     1005                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     1006 
     1007             IF(ll_tmp1) THEN 
     1008               zcpy(ji,jj) = 1.0_wp 
     1009             ELSE IF(ll_tmp2) THEN 
     1010               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
     1011               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
     1012                           &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     1013                zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 
     1014 
    9971015               ELSE 
    9981016                  zcpy(ji,jj) = 0._wp 
     
    11851203                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
    11861204               ENDIF 
    1187                IF( ln_wd ) THEN 
    1188                   zdpdx1 = zdpdx1 * zcpx(ji,jj) 
    1189                   zdpdx2 = zdpdx2 * zcpx(ji,jj) 
     1205               IF( ln_wd_il ) THEN 
     1206                  zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj) 
     1207                  zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 
    11901208               ENDIF 
    11911209               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
     
    12441262                  zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
    12451263               ENDIF 
    1246                IF( ln_wd ) THEN 
    1247                   zdpdy1 = zdpdy1 * zcpy(ji,jj) 
    1248                   zdpdy2 = zdpdy2 * zcpy(ji,jj) 
     1264               IF( ln_wd_il ) THEN 
     1265                  zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj)  
     1266                  zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj)  
    12491267               ENDIF 
    12501268 
     
    12561274      END DO 
    12571275      ! 
    1258       IF( ln_wd )   DEALLOCATE( zcpx, zcpy ) 
     1276      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy ) 
    12591277      ! 
    12601278   END SUBROUTINE hpg_prj 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r9019 r9023  
    2828   USE dom_oce        ! ocean space and time domain 
    2929   USE sbc_oce        ! Surface boundary condition: ocean fields 
     30   USE sbcrnf         ! river runoffs 
    3031   USE phycst         ! physical constants 
    3132   USE dynadv         ! dynamics: vector invariant versus flux form 
     
    130131            ! so that asselin contribution is removed at the same time  
    131132            DO jk = 1, jpkm1 
    132                un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 
    133                vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
     133               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk) 
     134               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
    134135            END DO   
    135136         ENDIF 
     
    207208            ! (used as a now filtered scale factor until the swap) 
    208209            ! ---------------------------------------------------- 
    209             IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN    ! No asselin filtering on thicknesses if forward time splitting 
    210                e3t_b(:,:,1:jpkm1) = e3t_n(:,:,1:jpkm1) 
    211             ELSE 
    212                DO jk = 1, jpkm1 
    213                   e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) 
     210            DO jk = 1, jpkm1 
     211               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) 
     212            END DO 
     213            ! Add volume filter correction: compatibility with tracer advection scheme 
     214            ! => time filter + conservation correction (only at the first level) 
     215            zcoef = atfp * rdt * r1_rau0 
     216            IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting 
     217               e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) & 
     218                              &                      - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 
     219            ELSE                     ! if ice shelf melting 
     220               DO jj = 1, jpj 
     221                  DO ji = 1, jpi 
     222                     ikt = mikt(ji,jj) 
     223                     e3t_b(ji,jj,ikt) = e3t_b(ji,jj,ikt) - zcoef * (  emp_b   (ji,jj) - emp   (ji,jj)  & 
     224                        &                                           - rnf_b   (ji,jj) + rnf   (ji,jj)  & 
     225                        &                                           + fwfisf_b(ji,jj) - fwfisf(ji,jj)  ) * tmask(ji,jj,ikt) 
     226                  END DO 
    214227               END DO 
    215228               ! Add volume filter correction: compatibility with tracer advection scheme 
     
    218231               IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting 
    219232                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) & 
    220                                  &                      - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 
     233                                 &                      ) * tmask(:,:,1) 
     234                  IF( ln_rnf_depth ) THEN  
     235                     DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too, not sure for ice shelf case yet 
     236                        DO jj = 1, jpj 
     237                           DO ji = 1, jpi 
     238                              IF( mikt(ji,jj) <= jk .and. jk <=  nk_rnf(ji,jj)  ) THEN 
     239                                 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef *  (-rnf_b(ji,jj) + rnf(ji,jj))*(e3t_n(ji,jj,jk)/h_rnf(ji,jj)  )*tmask(ji,jj,jk) 
     240                              ENDIF 
     241                           ENDDO 
     242                        ENDDO 
     243                     ENDDO 
     244                  ELSE 
     245                      e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 
     246                  ENDIF 
    221247               ELSE                     ! if ice shelf melting 
    222248                  DO jj = 1, jpj 
     
    229255                  END DO 
    230256               END IF 
    231             ENDIF 
     257            END IF 
    232258            ! 
    233259            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r9019 r9023  
    7474      ! 
    7575      INTEGER  ::   ji, jj, jk                   ! dummy loop indices 
    76       REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r   ! local scalars 
     76      REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r, zld   ! local scalars 
    7777      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zpice 
    7878      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
     
    121121               END DO  
    122122            END DO 
     123            ! 
     124            IF (ln_scal_load) THEN 
     125               zld = rn_scal_load * grav 
     126               DO jj = 2, jpjm1                    ! add scalar approximation for load potential 
     127                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     128                     spgu(ji,jj) = spgu(ji,jj) + zld * ( sshn(ji+1,jj) - sshn(ji,jj) ) * r1_e1u(ji,jj) 
     129                     spgv(ji,jj) = spgv(ji,jj) + zld * ( sshn(ji,jj+1) - sshn(ji,jj) ) * r1_e2v(ji,jj) 
     130                  END DO  
     131               END DO 
     132            ENDIF 
    123133         ENDIF 
    124134         ! 
     
    181191      NAMELIST/namdyn_spg/ ln_dynspg_exp       , ln_dynspg_ts,   & 
    182192      &                    ln_bt_fw, ln_bt_av  , ln_bt_auto  ,   & 
    183       &                    nn_baro , rn_bt_cmax, nn_bt_flt 
     193      &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha 
    184194      !!---------------------------------------------------------------------- 
    185195      ! 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r9019 r9023  
    11MODULE dynspg_ts 
     2 
     3   !! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !  
     4 
    25   !!====================================================================== 
    36   !!                   ***  MODULE  dynspg_ts  *** 
     
    5760   USE restart         ! only for lrst_oce 
    5861   USE timing          ! Timing     
     62   USE diatmb          ! Top,middle,bottom output 
     63#if defined key_agrif 
     64   USE agrif_opa_interp ! agrif 
     65   USE agrif_oce 
     66#endif 
     67#if defined key_asminc    
     68   USE asminc          ! Assimilation increment 
     69#endif 
    5970 
    6071   IMPLICIT NONE 
     
    6879   !! Time filtered arrays at baroclinic time step: 
    6980   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv   !: Advection vel. at "now" barocl. step 
    70  
    71    INTEGER , SAVE ::   icycle   ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro 
    72    REAL(wp), SAVE ::   rdtbt    ! Barotropic time step 
     81   INTEGER, SAVE :: icycle      ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro 
     82   REAL(wp),SAVE :: rdtbt       ! Barotropic time step 
    7383   ! 
    7484   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   wgtbtp1, wgtbtp2   ! 1st & 2nd weights used in time filtering of barotropic fields 
     
    131141      !!      -Update the filtered free surface at step "n+1"      : ssha 
    132142      !!      -Update filtered barotropic velocities at step "n+1" : ua_b, va_b 
    133       !!      -Compute barotropic advective velocities at step "n" : un_adv, vn_adv 
     143      !!      -Compute barotropic advective fluxes at step "n"    : un_adv, vn_adv 
    134144      !!      These are used to advect tracers and are compliant with discrete 
    135145      !!      continuity equation taken at the baroclinic time steps. This  
     
    159169      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points 
    160170      ! 
     171      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True.  
     172 
     173      INTEGER  :: iwdg, jwdg, kwdg   ! short-hand values for the indices of the output point 
     174 
     175      REAL(wp) ::   zepsilon, zgamma            !   -      - 
    161176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy   ! Wetting/Dying gravity filter coef. 
     177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 
     178      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2    ! averages over the sub-steps of zuwdmask and zvwdmask 
    162179      !!---------------------------------------------------------------------- 
    163180      ! 
    164181      IF( ln_timing )   CALL timing_start('dyn_spg_ts') 
    165182      ! 
    166       IF( ln_wd ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) ) 
     183      IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) ) 
     184      !                                         !* Allocate temporary arrays 
     185      IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj)) 
    167186      ! 
    168187      zmdi=1.e+20                               !  missing data indicator for masking 
    169188      ! 
    170       !                                            ! reciprocal of baroclinic time step  
     189      zwdramp = r_rn_wdmin1               ! simplest ramp  
     190!     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 
     191      !                                         ! reciprocal of baroclinic time step  
    171192      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt 
    172193      ELSE                                        ;   z2dt_bf = 2.0_wp * rdt 
     
    174195      z1_2dt_b = 1.0_wp / z2dt_bf  
    175196      ! 
    176       ll_init     = ln_bt_av                       ! if no time averaging, then no specific restart  
     197      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart  
    177198      ll_fw_start = .FALSE. 
    178       !                                            ! time offset in steps for bdy data update 
     199      !                                         ! time offset in steps for bdy data update 
    179200      IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_baro 
    180201      ELSE                       ;   noffset =   0  
    181202      ENDIF 
    182203      ! 
    183       IF( kt == nit000 ) THEN                !* initialisation 
     204      IF( kt == nit000 ) THEN                   !* initialisation 
    184205         ! 
    185206         IF(lwp) WRITE(numout,*) 
     
    405426      !                                   ! ---------------------------------------------------- 
    406427      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient 
    407          IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters 
    408             DO jj = 2, jpjm1 
    409                DO ji = 2, jpim1  
    410                   ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    411                      &      MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            & 
    412                      &      MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 
    413                      &                                                         > rn_wdmin1 + rn_wdmin2 
    414                   ll_tmp2 = ( ABS( sshn(ji+1,jj)             -   sshn(ji  ,jj))  > 1.E-12 ).AND.( & 
    415                      &      MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    416                      &      MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    417                      ! 
    418                   IF(ll_tmp1) THEN 
    419                      zcpx(ji,jj) = 1.0_wp 
    420                   ELSE IF(ll_tmp2) THEN   ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    421                      zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 
    422                         &             / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
    423                   ELSE 
    424                      zcpx(ji,jj) = 0._wp 
    425                   ENDIF 
    426                   ! 
    427                   ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    428                      &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            & 
    429                      &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 
    430                      &                                                         > rn_wdmin1 + rn_wdmin2 
    431                   ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1))  > 1.E-12 ).AND.( & 
    432                      &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    433                      &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    434                      ! 
    435                   IF(ll_tmp1) THEN 
    436                      zcpy(ji,jj) = 1.0_wp 
    437                   ELSE IF(ll_tmp2) THEN 
    438                      ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    439                      zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 
    440                         &             / (sshn(ji,jj+1) -  sshn(ji,jj  )) ) 
    441                   ELSE 
    442                      zcpy(ji,jj) = 0._wp 
    443                   ENDIF 
    444                END DO 
    445             END DO 
    446             ! 
    447             DO jj = 2, jpjm1 
    448                DO ji = 2, jpim1 
    449                   zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
    450                      &                            * r1_e1u(ji,jj) * zcpx(ji,jj) 
    451                   zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
    452                      &                            * r1_e2v(ji,jj) * zcpy(ji,jj) 
    453                END DO 
    454             END DO 
    455             ! 
     428        IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters 
     429           DO jj = 2, jpjm1 
     430              DO ji = 2, jpim1  
     431                ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                & 
     432                     &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
     433                     &    MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     434                     &                                                       > rn_wdmin1 + rn_wdmin2 
     435                ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( & 
     436                     &    MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     437                     &    MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     438    
     439                IF(ll_tmp1) THEN 
     440                  zcpx(ji,jj) = 1.0_wp 
     441                ELSE IF(ll_tmp2) THEN 
     442                  ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     443                  zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
     444                              &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
     445                  zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
     446 
     447                ELSE 
     448                  zcpx(ji,jj) = 0._wp 
     449                END IF 
     450          
     451                ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                & 
     452                     &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
     453                     &    MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
     454                     &                                                       > rn_wdmin1 + rn_wdmin2 
     455                ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( & 
     456                     &    MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     457                     &    MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     458   
     459                IF(ll_tmp1) THEN 
     460                  zcpy(ji,jj) = 1.0_wp 
     461                ELSE IF(ll_tmp2) THEN 
     462                  ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
     463                  zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
     464                              &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     465                  zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 
     466 
     467                ELSE 
     468                  zcpy(ji,jj) = 0._wp 
     469                END IF 
     470              END DO 
     471           END DO 
     472  
     473           DO jj = 2, jpjm1 
     474              DO ji = 2, jpim1 
     475                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
     476                        &                        * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth 
     477                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
     478                        &                        * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth 
     479 
     480              END DO 
     481           END DO 
     482 
    456483         ELSE 
    457484            ! 
     
    473500      END DO  
    474501      ! 
    475       !                 ! Add BOTTOM stress contribution from baroclinic velocities:       
    476       IF( ln_bt_fw ) THEN 
     502      !                                         ! Add bottom stress contribution from baroclinic velocities:       
     503      IF (ln_bt_fw) THEN 
    477504         DO jj = 2, jpjm1                           
    478505            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    495522      ! 
    496523      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
    497       IF( ln_wd ) THEN 
    498          zztmp = - 1._wp / rdtbt 
    499          DO jj = 2, jpjm1                           
     524      IF( ln_wd_il ) THEN 
     525        zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) *  wdrampu(ji,jj) 
     526        zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) *  wdrampv(ji,jj) 
     527      ELSE 
     528        zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 
     529        zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 
     530      END IF 
     531      ! 
     532      !                                         ! Add top stress contribution from baroclinic velocities:       
     533      IF( ln_bt_fw ) THEN 
     534         DO jj = 2, jpjm1 
    500535            DO ji = fs_2, fs_jpim1   ! vector opt. 
    501536               zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) * zwx(ji,jj) 
     
    657692      vn_adv(:,:) = 0._wp 
    658693      !                                             ! ==================== ! 
     694 
     695      IF (ln_wd_dl) THEN 
     696         zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary)  
     697         zvwdmask(:,:) = 0._wp  !  
     698         zuwdav2(:,:) =  0._wp  
     699         zvwdav2(:,:) =  0._wp    
     700      END IF  
     701 
     702 
    659703      DO jn = 1, icycle                             !  sub-time-step loop  ! 
    660704         !                                          ! ==================== ! 
     
    684728            ! Extrapolate Sea Level at step jit+0.5: 
    685729            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 
    686             ! 
     730             
     731            ! set wetting & drying mask at tracer points for this barotropic sub-step  
     732            IF ( ln_wd_dl ) THEN  
     733 
     734               IF ( ln_wd_dl_rmp ) THEN  
     735                  DO jj = 1, jpj                                  
     736                     DO ji = 1, jpi   ! vector opt.   
     737                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN  
     738!                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin2 ) THEN  
     739                           ztwdmask(ji,jj) = 1._wp 
     740                        ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN 
     741                           ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1)) )  
     742                        ELSE  
     743                           ztwdmask(ji,jj) = 0._wp 
     744                        END IF 
     745                     END DO 
     746                  END DO 
     747               ELSE 
     748                  DO jj = 1, jpj                                  
     749                     DO ji = 1, jpi   ! vector opt.   
     750                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN  
     751                           ztwdmask(ji,jj) = 1._wp 
     752                        ELSE  
     753                           ztwdmask(ji,jj) = 0._wp 
     754                        END IF 
     755                     END DO 
     756                  END DO 
     757               END IF  
     758 
     759            END IF  
     760            
     761 
    687762            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
    688763               DO ji = 2, fs_jpim1   ! Vector opt. 
     
    736811         ENDIF 
    737812#endif 
    738          IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    739          ! 
     813         IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
     814 
     815         IF ( ln_wd_dl ) THEN  
     816 
     817 
     818! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells  
     819 
     820            DO jj = 1, jpjm1                                  
     821               DO ji = 1, jpim1    
     822                  IF ( zwx(ji,jj) > 0.0 ) THEN  
     823                     zuwdmask(ji, jj) = ztwdmask(ji  ,jj)  
     824                  ELSE  
     825                     zuwdmask(ji, jj) = ztwdmask(ji+1,jj)   
     826                  END IF  
     827                  zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 
     828                  un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 
     829 
     830                  IF ( zwy(ji,jj) > 0.0 ) THEN 
     831                     zvwdmask(ji, jj) = ztwdmask(ji, jj  ) 
     832                  ELSE  
     833                     zvwdmask(ji, jj) = ztwdmask(ji, jj+1)   
     834                  END IF  
     835                  zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj)  
     836                  vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 
     837               END DO 
     838            END DO 
     839 
     840 
     841         END IF     
     842          
    740843         ! Sum over sub-time-steps to compute advective velocities 
    741844         za2 = wgtbtp2(jn) 
    742845         un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 
    743846         vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 
    744          ! 
     847          
     848         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)  
     849         IF ( ln_wd_dl_bc ) THEN 
     850            zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 
     851            zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 
     852         END IF  
     853 
    745854         ! Set next sea level: 
    746855         DO jj = 2, jpjm1                                  
     
    788897           za3= 0._wp               
    789898         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880  
    790            za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps     
    791            za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps  
    792            za2=0.088_wp                     ! za2 = gam 
    793            za3=0.013_wp                     ! za3 = eps 
     899            IF (rn_bt_alpha==0._wp) THEN 
     900               za0=0.614_wp                 ! za0 = 1/2 +   gam + 2*eps 
     901               za1=0.285_wp                 ! za1 = 1/2 - 2*gam - 3*eps 
     902               za2=0.088_wp                 ! za2 = gam 
     903               za3=0.013_wp                 ! za3 = eps 
     904            ELSE 
     905               zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 
     906               zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 
     907               za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 
     908               za1 = 1._wp - za0 - zgamma - zepsilon 
     909               za2 = zgamma 
     910               za3 = zepsilon 
     911            ENDIF  
    794912         ENDIF 
    795913         ! 
    796914         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) & 
    797915          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
    798          IF( ln_wd ) THEN                   ! Calculating and applying W/D gravity filters 
     916         IF( ln_wd_il ) THEN                   ! Calculating and applying W/D gravity filters 
    799917           DO jj = 2, jpjm1 
    800918              DO ji = 2, jpim1  
    801919                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                & 
    802                      &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji+1,jj) ) .AND.            & 
    803                      &    MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) ) & 
     920                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) .AND.            & 
     921                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 
    804922                     &                                                             > rn_wdmin1 + rn_wdmin2 
    805923                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji+1,jj))  > 1.E-12 ).AND.( & 
    806924                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                & 
    807                      &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     925                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    808926    
    809927                IF(ll_tmp1) THEN 
     
    811929                ELSE IF(ll_tmp2) THEN 
    812930                  ! no worries about  zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj) = 0, it won't happen ! here 
    813                   zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +    ht_wd(ji+1,jj) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) & 
     931                  zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +     ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 
    814932                              &    / (zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj)) ) 
    815933                ELSE 
     
    818936          
    819937                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
    820                      &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji,jj+1) ) .AND.            & 
    821                      &    MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) ) & 
     938                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) .AND.            & 
     939                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 
    822940                     &                                                             > rn_wdmin1 + rn_wdmin2 
    823941                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji,jj+1))  > 1.E-12 ).AND.( & 
    824942                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
    825                      &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     943                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    826944    
    827945                IF(ll_tmp1) THEN 
     
    829947                ELSE IF(ll_tmp2) THEN 
    830948                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here 
    831                   zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +    ht_wd(ji,jj+1) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) & 
     949                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +     ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 
    832950                              &    / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  )) ) 
    833951                ELSE 
     
    839957         ! 
    840958         ! Compute associated depths at U and V points: 
    841          IF( .NOT.ln_linssh  .AND. .NOT.ln_dynadv_vec ) THEN   !* Vector form 
     959         IF( .NOT.ln_linssh  .AND. .NOT.ln_dynadv_vec ) THEN     !* Vector form 
    842960            !                                         
    843961            DO jj = 2, jpjm1                             
     
    9151033         ENDIF 
    9161034         ! 
    917          DO jj = 2, jpjm1               ! Add top/bottom stresses: 
    918             DO ji = fs_2, fs_jpim1   ! vector opt. 
    919                zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 
    920                zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 
    921             END DO 
    922          END DO 
     1035         ! Add bottom stresses: 
     1036!jth do implicitly instead 
     1037         IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs 
     1038            zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 
     1039            zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 
     1040         ENDIF  
    9231041         ! 
    9241042         ! Surface pressure trend: 
    925          IF( ln_wd ) THEN 
     1043         IF( ln_wd_il ) THEN 
    9261044           DO jj = 2, jpjm1 
    9271045              DO ji = 2, jpim1  
     
    9291047                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 
    9301048                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 
    931                  zwx(ji,jj) = zu_spg * zcpx(ji,jj)  
    932                  zwy(ji,jj) = zv_spg * zcpy(ji,jj) 
     1049                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj)  
     1050                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj) 
    9331051              END DO 
    9341052           END DO 
     
    9391057                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 
    9401058                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 
    941                  zwx(ji,jj) = zu_spg 
    942                  zwy(ji,jj) = zv_spg 
     1059                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg 
     1060                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg 
    9431061              END DO 
    9441062           END DO 
     
    9471065         ! 
    9481066         ! Set next velocities: 
    949          IF( ln_dynadv_vec .OR. ln_linssh ) THEN   !* Vector form 
     1067         IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form 
    9501068            DO jj = 2, jpjm1 
    9511069               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    9611079                            &                                 + zv_frc(ji,jj) ) & 
    9621080                            &   ) * ssvmask(ji,jj) 
     1081  
     1082!jth implicit bottom friction: 
     1083                  IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 
     1084                     ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rdtbt * bfrua(ji,jj) * hur_e(ji,jj)) 
     1085                     va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rdtbt * bfrva(ji,jj) * hvr_e(ji,jj)) 
     1086                  ENDIF 
     1087 
    9631088               END DO 
    9641089            END DO 
    9651090            ! 
    966          ELSE                                      !* Flux form 
     1091         ELSE                           !* Flux form 
    9671092            DO jj = 2, jpjm1 
    9681093               DO ji = fs_2, fs_jpim1   ! vector opt. 
    9691094 
    970                   IF( ln_wd ) THEN 
    971                     zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 
    972                     zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 
    973                   ELSE 
    974                     zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
    975                     zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
    976                   END IF 
     1095                  zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
     1096                  zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
     1097 
    9771098                  zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 
    9781099                  zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 
     
    9921113            END DO 
    9931114         ENDIF 
    994          ! 
     1115 
     1116          
    9951117         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only) 
    996             IF( ln_wd ) THEN 
    997               hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 
    998               hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 
    999             ELSE 
    1000               hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
    1001               hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
    1002             END IF 
     1118            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     1119            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
    10031120            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 
    10041121            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 
     
    10331150            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:)  
    10341151            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:)  
    1035          ELSE                                              ! Sum transports 
    1036             ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
    1037             va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
    1038          ENDIF 
    1039          !                                   ! Sum sea level 
     1152         ELSE                                       ! Sum transports 
     1153            IF ( .NOT.ln_wd_dl ) THEN   
     1154               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
     1155               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
     1156            ELSE  
     1157               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:) 
     1158               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:) 
     1159            END IF  
     1160         ENDIF 
     1161         !                                          ! Sum sea level 
    10401162         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 
     1163 
    10411164         !                                                 ! ==================== ! 
    10421165      END DO                                               !        end loop      ! 
     
    10471170      ! 
    10481171      ! Set advection velocity correction: 
    1049       zwx(:,:) = un_adv(:,:) 
    1050       zwy(:,:) = vn_adv(:,:) 
    1051       IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN      
    1052          un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:) 
    1053          vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:) 
    1054       ELSE 
    1055          un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:) 
    1056          vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:) 
    1057       END IF 
    1058  
    1059       IF( ln_bt_fw ) THEN ! Save integrated transport for next computation 
     1172      IF (ln_bt_fw) THEN 
     1173         zwx(:,:) = un_adv(:,:) 
     1174         zwy(:,:) = vn_adv(:,:) 
     1175         IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 
     1176            un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) ) 
     1177            vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) ) 
     1178            ! 
     1179            ! Update corrective fluxes for next time step: 
     1180            un_bf(:,:)  = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:)) 
     1181            vn_bf(:,:)  = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:)) 
     1182         ELSE 
     1183            un_bf(:,:) = 0._wp 
     1184            vn_bf(:,:) = 0._wp  
     1185         END IF          
     1186         ! Save integrated transport for next computation 
    10601187         ub2_b(:,:) = zwx(:,:) 
    10611188         vb2_b(:,:) = zwy(:,:) 
    10621189      ENDIF 
     1190 
     1191 
    10631192      ! 
    10641193      ! Update barotropic trend: 
     
    10901219         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
    10911220      ENDIF 
    1092       ! 
     1221 
     1222 
     1223      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)   
    10931224      DO jk = 1, jpkm1 
    1094          ! Correct velocities: 
    1095          un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 
    1096          vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 
    1097          ! 
     1225         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk) 
     1226         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 
    10981227      END DO 
    1099       ! 
    1100       CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current 
    1101       CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic i-current 
     1228 
     1229      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN  
     1230         DO jk = 1, jpkm1 
     1231            un(:,:,jk) = ( un_adv(:,:) + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)) ) * umask(:,:,jk)  
     1232            vn(:,:,jk) = ( vn_adv(:,:) + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)) ) * vmask(:,:,jk)   
     1233         END DO 
     1234      END IF  
     1235 
     1236       
     1237      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu_n(:,:) )    ! barotropic i-current 
     1238      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv_n(:,:) )    ! barotropic i-current 
    11021239      ! 
    11031240#if defined key_agrif 
     
    11191256      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' ) 
    11201257      ! 
    1121       IF( ln_wd )   DEALLOCATE( zcpx, zcpy ) 
     1258      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy ) 
     1259      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 ) 
    11221260      ! 
    11231261      IF( ln_diatmb ) THEN 
     
    12221360         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )    
    12231361         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) )  
     1362         CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:) )    
     1363         CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:) )  
    12241364         IF( .NOT.ln_bt_av ) THEN 
    12251365            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )    
     
    12411381         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) ) 
    12421382         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) ) 
     1383         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:) ) 
     1384         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:) ) 
    12431385         ! 
    12441386         IF (.NOT.ln_bt_av) THEN 
     
    12911433      rdtbt = rdt / REAL( nn_baro , wp ) 
    12921434      zcmax = zcmax * rdtbt 
    1293                      ! Print results 
     1435      ! Print results 
    12941436      IF(lwp) WRITE(numout,*) 
    12951437      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 
     
    13171459#if defined key_agrif 
    13181460      ! Restrict the use of Agrif to the forward case only 
    1319       IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 
     1461!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 
    13201462#endif 
    13211463      ! 
     
    13331475      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax 
    13341476      ! 
     1477      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha 
     1478      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN 
     1479         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' ) 
     1480      ENDIF 
     1481      ! 
    13351482      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN 
    13361483         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' ) 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r9019 r9023  
    2727   USE agrif_opa_interp 
    2828#endif 
    29 #if defined key_asminc    
    30    USE   asminc       ! Assimilation increment 
    31 #endif 
    3229   ! 
    3330   USE in_out_manager ! I/O manager 
     
    3734   USE lib_mpp        ! MPP library 
    3835   USE timing         ! Timing 
    39    USE wet_dry        ! Wetting/Drying flux limting 
     36   USE wet_dry        ! Wetting/Drying flux limiting 
    4037 
    4138   IMPLICIT NONE 
     
    9188      !                                           !   After Sea Surface Height   ! 
    9289      !                                           !------------------------------! 
    93       IF(ln_wd) THEN 
     90      IF(ln_wd_il) THEN 
    9491         CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 
    9592      ENDIF 
     
    106103      !  
    107104      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
    108  
     105      ! 
     106#if defined key_agrif 
     107      CALL agrif_ssh( kt ) 
     108#endif 
     109      ! 
    109110      IF ( .NOT.ln_dynspg_ts ) THEN 
    110          ! These lines are not necessary with time splitting since 
    111          ! boundary condition on sea level is set during ts loop 
    112 # if defined key_agrif 
    113          CALL agrif_ssh( kt ) 
    114 # endif 
    115111         IF( ln_bdy ) THEN 
    116112            CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary 
     
    118114         ENDIF 
    119115      ENDIF 
    120  
    121 #if defined key_asminc 
    122       IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN     ! Include the IAU weighted SSH increment 
    123          CALL ssh_asm_inc( kt ) 
    124          ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:) 
    125       ENDIF 
    126 #endif 
    127116      !                                           !------------------------------! 
    128117      !                                           !           outputs            ! 
     
    209198      ENDIF 
    210199      ! 
    211       IF( ln_timing )   CALL timing_stop('wzv') 
     200#if defined key_agrif  
     201      IF( .NOT. AGRIF_Root() ) THEN  
     202         IF ((nbondi ==  1).OR.(nbondi == 2)) wn(nlci-1 , :     ,:) = 0.e0      ! east  
     203         IF ((nbondi == -1).OR.(nbondi == 2)) wn(2      , :     ,:) = 0.e0      ! west  
     204         IF ((nbondj ==  1).OR.(nbondj == 2)) wn(:      ,nlcj-1 ,:) = 0.e0      ! north  
     205         IF ((nbondj == -1).OR.(nbondj == 2)) wn(:      ,2      ,:) = 0.e0      ! south  
     206      ENDIF  
     207#endif  
     208      ! 
     209      IF( nn_timing == 1 )  CALL timing_stop('wzv') 
    212210      ! 
    213211   END SUBROUTINE wzv 
     
    246244      ENDIF 
    247245      !              !==  Euler time-stepping: no filter, just swap  ==! 
    248       IF(  ( neuler == 0 .AND. kt == nit000 ) .OR.    & 
    249          & ( ln_bt_fw    .AND. ln_dynspg_ts )      ) THEN  
     246      IF  ( neuler == 0 .AND. kt == nit000 ) THEN 
    250247         sshb(:,:) = sshn(:,:)                              ! before <-- now 
    251248         sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now) 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90

    r9019 r9023  
    11MODULE wet_dry 
     2 
     3   !! includes updates to namelist namwad for diagnostic outputs of ROMS wetting and drying 
     4 
    25   !!============================================================================== 
    36   !!                       ***  MODULE  wet_dry  *** 
    47   !! Wetting and drying includes initialisation routine and routines to 
    58   !! compute and apply flux limiters and preserve water depth positivity 
    6    !! only effects if wetting/drying is on (ln_wd == .true.) 
     9   !! only effects if wetting/drying is on (ln_wd_il == .true. or ln_wd_dl==.true. ) 
    710   !!============================================================================== 
    811   !! History :  3.6  ! 2014-09  ((H.Liu)  Original code 
     
    3336 
    3437   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask   !: u- and v- limiter  
    35    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   ht_wd    !: wetting and drying t-pt depths 
    3638   !                                                           !  (can include negative depths) 
    37  
    38    LOGICAL,  PUBLIC  ::   ln_wd       !: Wetting/drying activation switch (T:on,F:off) 
     39   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdramp, wdrampu, wdrampv !: for hpg limiting 
     40 
     41   LOGICAL,  PUBLIC  ::   ln_wd_il    !: Wetting/drying il activation switch (T:on,F:off) 
     42   LOGICAL,  PUBLIC  ::   ln_wd_dl    !: Wetting/drying dl activation switch (T:on,F:off) 
     43   REAL(wp), PUBLIC  ::   rn_wdmin0   !: depth at which wetting/drying starts 
    3944   REAL(wp), PUBLIC  ::   rn_wdmin1   !: minimum water depth on dried cells 
    40    REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells 
     45   REAL(wp), PUBLIC  ::   r_rn_wdmin1 !: 1/minimum water depth on dried cells  
     46   REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerance of minimum water depth on dried cells 
    4147   REAL(wp), PUBLIC  ::   rn_wdld     !: land elevation below which wetting/drying will be considered 
    4248   INTEGER , PUBLIC  ::   nn_wdit     !: maximum number of iteration for W/D limiter 
     49   LOGICAL,  PUBLIC  ::   ln_wd_dl_bc !: DL scheme: True implies 3D velocities are set to the barotropic values at points  
     50                                      !: where the flow is from wet points on less than half the barotropic sub-steps   
     51   LOGICAL,  PUBLIC  ::  ln_wd_dl_rmp !: use a ramp for the dl flux limiter between 2 rn_wdmin1 and rn_wdmin1 (rather than a cut-off at rn_wdmin1)       
     52   REAL(wp), PUBLIC  ::   ssh_ref     !: height of z=0 with respect to the geoid;  
     53 
     54   LOGICAL,  PUBLIC  ::   ll_wd       !: Wetting/drying activation switch if either ln_wd_il or ln_wd_dl 
    4355 
    4456   PUBLIC   wad_init                  ! initialisation routine called by step.F90 
     
    5971      !! ** input   : - namwad namelist 
    6072      !!---------------------------------------------------------------------- 
    61       INTEGER  ::   ios, ierr   ! Local integer 
    62       !! 
    63       NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 
     73      !! 
     74      NAMELIST/namwad/ ln_wd_il, ln_wd_dl, rn_wdmin0, rn_wdmin1, rn_wdmin2, rn_wdld, & 
     75                     & nn_wdit, ln_wd_dl_bc, ln_wd_dl_rmp 
     76      INTEGER  ::   ios                 ! Local integer output status for namelist read 
     77      INTEGER  ::   ierr                ! Local integer status array allocation  
    6478      !!---------------------------------------------------------------------- 
    6579      ! 
     
    7791         WRITE(numout,*) '~~~~~~~~' 
    7892         WRITE(numout,*) '   Namelist namwad' 
    79          WRITE(numout,*) '      Logical activation                 ln_wd      = ', ln_wd 
     93         WRITE(numout,*) '      Logical for Iter Lim wd option   ln_wd_il     = ', ln_wd_il 
     94         WRITE(numout,*) '      Logical for Dir. Lim wd option   ln_wd_dl     = ', ln_wd_dl 
     95         WRITE(numout,*) '      Depth at which wet/drying starts rn_wdmin0    = ', rn_wdmin0 
    8096         WRITE(numout,*) '      Minimum wet depth on dried cells rn_wdmin1    = ', rn_wdmin1 
    81          WRITE(numout,*) '      Tolerance of min wet depth     rn_wdmin2      = ', rn_wdmin2 
     97         WRITE(numout,*) '      Tolerance of min wet depth       rn_wdmin2    = ', rn_wdmin2 
    8298         WRITE(numout,*) '      land elevation threshold         rn_wdld      = ', rn_wdld 
    8399         WRITE(numout,*) '      Max iteration for W/D limiter    nn_wdit      = ', nn_wdit 
     100         WRITE(numout,*) '      T => baroclinic u,v=0 at dry pts: ln_wd_dl_bc = ', ln_wd_dl_bc      
     101         WRITE(numout,*) '      use a ramp for rwd limiter:  ln_wd_dl_rwd_rmp = ', ln_wd_dl_rmp 
     102  
    84103      ENDIF 
    85       ! 
    86       IF(ln_wd) THEN 
    87          ALLOCATE( wdmask(jpi,jpj), ht_wd(jpi,jpj),  STAT=ierr ) 
     104      IF( .NOT. ln_read_cfg ) THEN 
     105         WRITE(numout,*) '      No configuration file so seting ssh_ref to zero  ' 
     106          ssh_ref=0.0 
     107      ENDIF 
     108 
     109      r_rn_wdmin1=1/rn_wdmin1 
     110      ll_wd = .FALSE. 
     111      IF(ln_wd_il .OR. ln_wd_dl) THEN 
     112         ll_wd = .TRUE. 
     113         ALLOCATE( wdmask(jpi,jpj),   STAT=ierr ) 
     114         ALLOCATE( wdramp(jpi,jpj), wdrampu(jpi,jpj), wdrampv(jpi,jpj), STAT=ierr )  
    88115         IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 
    89116      ENDIF 
     
    118145      !!---------------------------------------------------------------------- 
    119146      ! 
    120       IF( ln_timing )   CALL timing_start('wad_lmt') 
    121       ! 
    122       !IF(lwp) WRITE(numout,*) 
    123       !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting' 
    124       ! 
     147      IF( nn_timing == 1 )  CALL timing_start('wad_lmt') 
     148      ! 
     149        
     150      DO jk = 1, jpkm1 
     151         un(:,:,jk) = un(:,:,jk) * zwdlmtu(:,:)  
     152         vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:,:)  
     153      END DO 
    125154      jflag  = 0 
    126155      zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen 
    127       !  
     156 
    128157      zflxp(:,:)   = 0._wp 
    129158      zflxn(:,:)   = 0._wp 
    130159      zflxu(:,:)   = 0._wp 
    131160      zflxv(:,:)   = 0._wp 
    132       ! 
    133       zwdlmtu(:,:) = 1._wp 
    134       zwdlmtv(:,:) = 1._wp 
    135       !  
     161 
     162      zwdlmtu(:,:)  = 1._wp 
     163      zwdlmtv(:,:)  = 1._wp 
     164       
    136165      ! Horizontal Flux in u and v direction 
    137166      DO jk = 1, jpkm1   
     
    143172         END DO   
    144173      END DO 
    145       ! 
     174        
    146175      zflxu(:,:) = zflxu(:,:) * e2u(:,:) 
    147176      zflxv(:,:) = zflxv(:,:) * e1v(:,:) 
    148       !  
     177       
    149178      wdmask(:,:) = 1 
    150179      DO jj = 2, jpj 
    151180         DO ji = 2, jpi  
    152             ! 
    153             IF( tmask(ji, jj, 1) < 0.5_wp )   CYCLE   ! we don't care about land cells 
    154             IF( ht_wd(ji,jj)     > zdepwd )  CYCLE   ! and cells which are unlikely to dry 
    155             ! 
    156             zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj) , 0._wp )  & 
    157                &         + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji,jj-1) , 0._wp )  
    158             zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj) , 0._wp )  & 
    159                &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,jj-1) , 0._wp )  
    160             ! 
    161             zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 
    162             IF( zdep2 .le. 0._wp) THEN  !add more safty, but not necessary 
    163                sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 
    164                IF( zflxu(ji,  jj) > 0._wp )  zwdlmtu(ji  ,jj) = 0._wp 
    165                IF( zflxu(ji-1,jj) < 0._wp )  zwdlmtu(ji-1,jj) = 0._wp 
    166                IF( zflxv(ji,  jj) > 0._wp )  zwdlmtv(ji  ,jj) = 0._wp 
    167                IF( zflxv(ji,jj-1) < 0._wp )  zwdlmtv(ji,jj-1) = 0._wp  
     181 
     182            IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE              ! we don't care about land cells 
     183            IF( ht_0(ji,jj) - ssh_ref > zdepwd ) CYCLE   ! and cells which are unlikely to dry 
     184 
     185            zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + & 
     186                         & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp)  
     187            zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + & 
     188                         & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp)  
     189 
     190            zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 
     191            IF(zdep2 .le. 0._wp) THEN  !add more safty, but not necessary 
     192               sshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 
     193               IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp 
     194               IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 
     195               IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp 
     196               IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp  
    168197               wdmask(ji,jj) = 0._wp 
    169             ENDIF 
     198            END IF 
     199         ENDDO 
     200      END DO 
     201 
     202 
     203! HPG limiter from jholt 
     204      wdramp(:,:) = min((ht_0(:,:) + sshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp) 
     205!jth assume don't need a lbc_lnk here 
     206      DO jj = 1, jpjm1 
     207         DO ji = 1, jpim1 
     208            wdrampu(ji,jj) = min(wdramp(ji,jj),wdramp(ji+1,jj)) 
     209            wdrampv(ji,jj) = min(wdramp(ji,jj),wdramp(ji,jj+1)) 
    170210         END DO 
    171211      END DO 
    172       !! 
    173       !! start limiter iterations  
     212! end HPG limiter 
     213 
     214 
     215       
     216        !! start limiter iterations  
    174217      DO jk1 = 1, nn_wdit + 1 
    175          !  
     218        
     219           
    176220         zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 
    177221         zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 
    178222         jflag = 0     ! flag indicating if any further iterations are needed 
    179          !  
     223           
    180224         DO jj = 2, jpj 
    181225            DO ji = 2, jpi  
    182                ! 
    183                IF( tmask(ji,jj,1) < 0.5_wp )  CYCLE  
    184                IF( ht_wd(ji,jj)   > zdepwd )   CYCLE 
    185                ! 
     226         
     227               IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE  
     228               IF( ht_0(ji,jj) > zdepwd )      CYCLE 
     229         
    186230               ztmp = e1e2t(ji,jj) 
    187                ! 
    188                zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj) , 0._wp )  & 
    189                   &   + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji,jj-1) , 0._wp )  
    190                zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj) , 0._wp )  & 
    191                   &   + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji,jj-1) , 0._wp )  
    192                ! 
     231 
     232               zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + & 
     233                      & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp)  
     234               zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + & 
     235                      & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp)  
     236           
    193237               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 
    194                zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 
    195                ! 
     238               zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 
     239           
    196240               IF( zdep1 > zdep2 ) THEN 
    197241                  wdmask(ji, jj) = 0 
     
    201245                  ! changes have zeroed the coefficient since further iterations will 
    202246                  ! not change anything 
    203                   IF( zcoef > 0._wp ) THEN   ;   jflag = 1  
    204                   ELSE                       ;   zcoef = 0._wp 
     247                  IF( zcoef > 0._wp ) THEN 
     248                     jflag = 1  
     249                  ELSE 
     250                     zcoef = 0._wp 
    205251                  ENDIF 
    206                   IF( jk1 > nn_wdit )  zcoef = 0._wp 
    207                   IF( zflxu1(ji,  jj) > 0._wp )  zwdlmtu(ji  ,jj) = zcoef 
    208                   IF( zflxu1(ji-1,jj) < 0._wp )  zwdlmtu(ji-1,jj) = zcoef 
    209                   IF( zflxv1(ji,  jj) > 0._wp )  zwdlmtv(ji  ,jj) = zcoef 
    210                   IF( zflxv1(ji,jj-1) < 0._wp )  zwdlmtv(ji,jj-1) = zcoef 
    211                ENDIF 
     252                  IF(jk1 > nn_wdit) zcoef = 0._wp 
     253                  IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef 
     254                  IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 
     255                  IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef 
     256                  IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 
     257               END IF 
    212258            END DO ! ji loop 
    213259         END DO  ! jj loop 
    214          ! 
     260 
    215261         CALL lbc_lnk( zwdlmtu, 'U', 1. ) 
    216262         CALL lbc_lnk( zwdlmtv, 'V', 1. ) 
    217          ! 
    218          IF(lk_mpp)   CALL mpp_max(jflag)   !max over the global domain 
    219          ! 
    220          IF(jflag == 0)   EXIT 
    221          !  
     263 
     264         IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain 
     265 
     266         IF(jflag == 0) EXIT 
     267           
    222268      END DO  ! jk1 loop 
    223269        
    224270      DO jk = 1, jpkm1 
    225          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:,:)  
    226          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:,:)  
    227       END DO 
    228  
    229 !!gm ==> Andrew  : the lbclnk below is useless since above lbclnk is applied on zwdlmtu/v 
    230 !!                                             and un, vn always with lbclnk 
     271         un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :)  
     272         vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :)  
     273      END DO 
     274 
    231275      CALL lbc_lnk( un, 'U', -1. ) 
    232276      CALL lbc_lnk( vn, 'V', -1. ) 
    233 !!gm end 
    234       ! 
    235       un_b(:,:) = un_b(:,:) * zwdlmtu(:,:) 
    236       vn_b(:,:) = vn_b(:,:) * zwdlmtv(:,:) 
    237 !!gm ==> Andrew   : probably same as above 
     277        ! 
     278      un_b(:,:) = un_b(:,:) * zwdlmtu(:, :) 
     279      vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :) 
    238280      CALL lbc_lnk( un_b, 'U', -1. ) 
    239281      CALL lbc_lnk( vn_b, 'V', -1. ) 
    240 !!gm end 
    241282        
    242283      IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 
     
    245286      !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field) 
    246287      ! 
     288      ! 
     289      ! 
     290      ! 
     291      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt') 
    247292      ! 
    248293      IF( ln_timing )   CALL timing_stop('wad_lmt') 
     
    276321      !!---------------------------------------------------------------------- 
    277322      ! 
    278       IF( ln_timing )  CALL timing_start('wad_lmt_bt') 
    279       !        
    280       !IF(lwp) WRITE(numout,*) 
    281       !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting' 
    282         
     323      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt') 
     324      ! 
    283325      jflag  = 0 
    284326      zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes 
     
    296338      DO jj = 2, jpj 
    297339         DO ji = 2, jpi  
    298             ! 
    299             IF( tmask(ji,jj,1) < 0.5_wp )  CYCLE   ! we don't care about land cells 
    300             IF( ht_wd(ji,jj)   > zdepwd )   CYCLE   ! and cells which are unlikely to dry 
    301             ! 
    302             zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj) , 0._wp )  & 
    303                &         + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji,jj-1) , 0._wp )  
    304             zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj) , 0._wp )  & 
    305                &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,jj-1) , 0._wp )  
    306  
    307             zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 
     340 
     341           IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE   ! we don't care about land cells 
     342           IF( ht_0(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry 
     343 
     344            zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + & 
     345                         & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp)  
     346            zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + & 
     347                         & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp)  
     348 
     349            zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 
    308350            IF(zdep2 .le. 0._wp) THEN  !add more safety, but not necessary 
    309                sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 
    310                IF( zflxu(ji,  jj) > 0._wp )  zwdlmtu(ji  ,jj) = 0._wp 
    311                IF( zflxu(ji-1,jj) < 0._wp )  zwdlmtu(ji-1,jj) = 0._wp 
    312                IF( zflxv(ji,  jj) > 0._wp )  zwdlmtv(ji  ,jj) = 0._wp 
    313                IF( zflxv(ji,jj-1) < 0._wp )  zwdlmtv(ji,jj-1) = 0._wp  
    314             ENDIF 
    315          END DO 
     351              sshn_e(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 
     352              IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp 
     353              IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 
     354              IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp 
     355              IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp  
     356            END IF 
     357         ENDDO 
    316358      END DO 
    317359 
     
    319361      !! start limiter iterations  
    320362      DO jk1 = 1, nn_wdit + 1 
    321          !  
     363        
     364           
    322365         zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 
    323366         zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 
    324367         jflag = 0     ! flag indicating if any further iterations are needed 
    325          ! 
     368           
    326369         DO jj = 2, jpj 
    327370            DO ji = 2, jpi  
    328                ! 
    329                IF( tmask(ji,jj, 1 ) < 0.5_wp  )  CYCLE  
    330                IF( ht_wd(ji,jj)      > zdepwd )   CYCLE 
    331                ! 
     371         
     372               IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE  
     373               IF( ht_0(ji,jj) > zdepwd )      CYCLE 
     374         
    332375               ztmp = e1e2t(ji,jj) 
    333                ! 
    334                zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj) , 0._wp )  & 
    335                   &   + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji,jj-1) , 0._wp )  
    336                zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj) , 0._wp )  & 
    337                   &   + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji,jj-1) , 0._wp )  
     376 
     377               zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + & 
     378                      & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp)  
     379               zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + & 
     380                      & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp)  
    338381           
    339382               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 
    340                zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 
     383               zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 
    341384           
    342385               IF(zdep1 > zdep2) THEN 
    343                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
    344                   !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 
    345                   ! flag if the limiter has been used but stop flagging if the only 
    346                   ! changes have zeroed the coefficient since further iterations will 
    347                   ! not change anything 
    348                   IF( zcoef > 0._wp ) THEN 
    349                       jflag = 1  
    350                   ELSE 
    351                       zcoef = 0._wp 
    352                   ENDIF 
    353                   IF( jk1 > nn_wdit )  zcoef = 0._wp 
    354                   IF( zflxu1(ji,  jj) > 0._wp )  zwdlmtu(ji  ,jj) = zcoef 
    355                   IF( zflxu1(ji-1,jj) < 0._wp )  zwdlmtu(ji-1,jj) = zcoef 
    356                   IF( zflxv1(ji,  jj) > 0._wp )  zwdlmtv(ji  ,jj) = zcoef 
    357                   IF( zflxv1(ji,jj-1) < 0._wp )  zwdlmtv(ji,jj-1) = zcoef 
    358                ENDIF 
     386                 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
     387                 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 
     388                 ! flag if the limiter has been used but stop flagging if the only 
     389                 ! changes have zeroed the coefficient since further iterations will 
     390                 ! not change anything 
     391                 IF( zcoef > 0._wp ) THEN 
     392                    jflag = 1  
     393                 ELSE 
     394                    zcoef = 0._wp 
     395                 ENDIF 
     396                 IF(jk1 > nn_wdit) zcoef = 0._wp 
     397                 IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef 
     398                 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 
     399                 IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef 
     400                 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 
     401               END IF 
    359402            END DO ! ji loop 
    360403         END DO  ! jj loop 
    361          ! 
     404 
    362405         CALL lbc_lnk( zwdlmtu, 'U', 1. ) 
    363406         CALL lbc_lnk( zwdlmtv, 'V', 1. ) 
    364          ! 
    365          IF(lk_mpp)   CALL mpp_max(jflag)   !max over the global domain 
    366          ! 
    367          IF( jflag == 0 )  EXIT 
    368          !     
     407 
     408         IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain 
     409 
     410         IF(jflag == 0) EXIT 
     411           
    369412      END DO  ! jk1 loop 
    370       !  
     413       
    371414      zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :)  
    372415      zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :)  
    373       ! 
     416 
    374417      CALL lbc_lnk( zflxu, 'U', -1. ) 
    375418      CALL lbc_lnk( zflxv, 'V', -1. ) 
    376       ! 
     419        
    377420      IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 
    378421        
     
    380423      !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field) 
    381424      ! 
    382       IF( ln_timing )  CALL timing_stop('wad_lmt_bt') 
    383       ! 
     425      ! 
     426      CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 
     427      CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 
     428      ! 
     429 
     430      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt') 
    384431   END SUBROUTINE wad_lmt_bt 
    385432 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r9019 r9023  
    4646   LOGICAL, PUBLIC :: ln_diaobs   !: Logical switch for the obs operator 
    4747   LOGICAL :: ln_sstnight         !: Logical switch for night mean SST obs 
    48     
     48   LOGICAL :: ln_sla_fp_indegs    !: T=> SLA obs footprint size specified in degrees, F=> in metres 
     49   LOGICAL :: ln_sst_fp_indegs    !: T=> SST obs footprint size specified in degrees, F=> in metres 
     50   LOGICAL :: ln_sss_fp_indegs    !: T=> SSS obs footprint size specified in degrees, F=> in metres 
     51   LOGICAL :: ln_sic_fp_indegs    !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 
     52 
     53   REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres) 
     54   REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres) 
     55   REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres) 
     56   REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres) 
     57   REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres) 
     58   REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres) 
     59   REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres) 
     60   REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres) 
     61 
    4962   INTEGER :: nn_1dint       !: Vertical interpolation method 
    50    INTEGER :: nn_2dint       !: Horizontal interpolation method 
     63   INTEGER :: nn_2dint       !: Default horizontal interpolation method 
     64   INTEGER :: nn_2dint_sla   !: SLA horizontal interpolation method  
     65   INTEGER :: nn_2dint_sst   !: SST horizontal interpolation method  
     66   INTEGER :: nn_2dint_sss   !: SSS horizontal interpolation method  
     67   INTEGER :: nn_2dint_sic   !: Seaice horizontal interpolation method  
    5168   INTEGER, DIMENSION(imaxavtypes) ::   nn_profdavtypes   !: Profile data types representing a daily average 
    5269   INTEGER :: nproftypes     !: Number of profile obs types 
    5370   INTEGER :: nsurftypes     !: Number of surface obs types 
    54    INTEGER, DIMENSION(:), ALLOCATABLE ::   nvarsprof, nvarssurf   !: Number of profile & surface variables 
    55    INTEGER, DIMENSION(:), ALLOCATABLE ::   nextrprof, nextrsurf   !: Number of profile & surface extra variables 
    56    INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sstbias_type   !: SST bias type     
    57    TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) ::   surfdata, surfdataqc   !: Initial surface data before & after quality control 
    58    TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) ::   profdata, profdataqc   !: Initial profile data before & after quality control 
     71   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     72      & nvarsprof, &         !: Number of profile variables 
     73      & nvarssurf            !: Number of surface variables 
     74   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     75      & nextrprof, &         !: Number of profile extra variables 
     76      & nextrsurf            !: Number of surface extra variables 
     77   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     78      & n2dintsurf           !: Interpolation option for surface variables 
     79   REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
     80      & zavglamscl, &        !: E/W diameter of averaging footprint for surface variables 
     81      & zavgphiscl           !: N/S diameter of averaging footprint for surface variables 
     82   LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
     83      & lfpindegs, &         !: T=> surface obs footprint size specified in degrees, F=> in metres 
     84      & llnightav            !: Logical for calculating night-time averages 
     85 
     86   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 
     87      & surfdata, &          !: Initial surface data 
     88      & surfdataqc           !: Surface data after quality control 
     89   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & 
     90      & profdata, &          !: Initial profile data 
     91      & profdataqc           !: Profile data after quality control 
    5992 
    6093   CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE ::   cobstypesprof, cobstypessurf   !: Profile & surface obs types 
     
    95128         & cn_profbfiles, &      ! T/S profile input filenames 
    96129         & cn_sstfbfiles, &      ! Sea surface temperature input filenames 
     130         & cn_sssfbfiles, &      ! Sea surface salinity input filenames 
    97131         & cn_slafbfiles, &      ! Sea level anomaly input filenames 
    98132         & cn_sicfbfiles, &      ! Seaice concentration input filenames 
    99133         & cn_velfbfiles, &      ! Velocity profile input filenames 
    100          & cn_sstbias_files      ! SST bias input filenames 
     134         & cn_sstbiasfiles      ! SST bias input filenames 
    101135      CHARACTER(LEN=128) :: & 
    102136         & cn_altbiasfile        ! Altimeter bias input filename 
     
    109143      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies  
    110144      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature 
     145      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity 
    111146      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration 
    112147      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
     
    116151      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files 
    117152      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
     153      LOGICAL :: ln_bound_reject ! Logical to remove obs near boundaries in LAMs. 
    118154      LOGICAL :: llvar1          ! Logical for profile variable 1 
    119155      LOGICAL :: llvar2          ! Logical for profile variable 1 
    120       LOGICAL :: llnightav       ! Logical for calculating night-time averages 
    121156      LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files 
    122157 
     
    134169 
    135170      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
    136          &            ln_sst, ln_sic, ln_vel3d,                       & 
    137          &            ln_altbias, ln_nea, ln_grid_global,             & 
    138          &            ln_grid_search_lookup,                          & 
    139          &            ln_ignmis, ln_s_at_t, ln_sstnight,              & 
     171         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               & 
     172         &            ln_altbias, ln_sstbias, ln_nea,                 & 
     173         &            ln_grid_global, ln_grid_search_lookup,          & 
     174         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
     175         &            ln_sstnight,                                    & 
     176         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
     177         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
    140178         &            cn_profbfiles, cn_slafbfiles,                   & 
    141179         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
    142          &            cn_velfbfiles, cn_altbiasfile,                  & 
     180         &            cn_velfbfiles, cn_sssfbfiles,                   & 
     181         &            cn_sstbiasfiles, cn_altbiasfile,                & 
    143182         &            cn_gridsearchfile, rn_gridsearchres,            & 
    144          &            rn_dobsini, rn_dobsend, nn_1dint, nn_2dint,     & 
     183         &            rn_dobsini, rn_dobsend,                         & 
     184         &            rn_sla_avglamscl, rn_sla_avgphiscl,             & 
     185         &            rn_sst_avglamscl, rn_sst_avgphiscl,             & 
     186         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
     187         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
     188         &            nn_1dint, nn_2dint,                             & 
     189         &            nn_2dint_sla, nn_2dint_sst,                     & 
     190         &            nn_2dint_sss, nn_2dint_sic,                     & 
    145191         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
    146          &            nn_profdavtypes, ln_sstbias, cn_sstbias_files 
     192         &            nn_profdavtypes 
    147193 
    148194      INTEGER :: jnumsstbias 
     
    157203      ! Read namelist parameters 
    158204      !----------------------------------------------------------------------- 
    159        
    160       !Initalise all values in namelist arrays 
    161       ALLOCATE(sstbias_type(jpmaxnfiles)) 
    162205      ! Some namelist arrays need initialising 
    163206      cn_profbfiles(:) = '' 
     
    166209      cn_sicfbfiles(:) = '' 
    167210      cn_velfbfiles(:) = '' 
    168       cn_sstbias_files(:) = '' 
     211      cn_sssfbfiles(:)    = '' 
     212      cn_sstbiasfiles(:) = '' 
    169213      nn_profdavtypes(:) = -1 
    170214 
     
    187231         RETURN 
    188232      ENDIF 
    189        
    190       !----------------------------------------------------------------------- 
    191       ! Set up list of observation types to be used 
    192       ! and the files associated with each type 
    193       !----------------------------------------------------------------------- 
    194  
    195       nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
    196       nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) ) 
    197  
    198       IF (ln_sstbias) THEN  
    199          lmask(:) = .FALSE.  
    200          WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE.  
    201          jnumsstbias = COUNT(lmask)  
    202          lmask(:) = .FALSE.  
    203       ENDIF       
    204  
    205       IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
    206          IF(lwp) WRITE(numout,cform_war) 
    207          IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
    208             &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 
    209             &                    ' are set to .FALSE. so turning off calls to dia_obs' 
    210          nwarn = nwarn + 1 
    211          ln_diaobs = .FALSE. 
    212          RETURN 
    213       ENDIF 
    214  
    215       IF ( nproftypes > 0 ) THEN 
    216  
    217          ALLOCATE( cobstypesprof(nproftypes) ) 
    218          ALLOCATE( ifilesprof(nproftypes) ) 
    219          ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
    220  
    221          jtype = 0 
    222          IF (ln_t3d .OR. ln_s3d) THEN 
    223             jtype = jtype + 1 
    224             clproffiles(jtype,:) = cn_profbfiles(:) 
    225             cobstypesprof(jtype) = 'prof  ' 
    226             ifilesprof(jtype) = 0 
    227             DO jfile = 1, jpmaxnfiles 
    228                IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 
    229                   ifilesprof(jtype) = ifilesprof(jtype) + 1 
    230             END DO 
    231          ENDIF 
    232          IF (ln_vel3d) THEN 
    233             jtype = jtype + 1 
    234             clproffiles(jtype,:) = cn_velfbfiles(:) 
    235             cobstypesprof(jtype) = 'vel   ' 
    236             ifilesprof(jtype) = 0 
    237             DO jfile = 1, jpmaxnfiles 
    238                IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 
    239                   ifilesprof(jtype) = ifilesprof(jtype) + 1 
    240             END DO 
    241          ENDIF 
    242  
    243       ENDIF 
    244  
    245       IF ( nsurftypes > 0 ) THEN 
    246  
    247          ALLOCATE( cobstypessurf(nsurftypes) ) 
    248          ALLOCATE( ifilessurf(nsurftypes) ) 
    249          ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
    250  
    251          jtype = 0 
    252          IF (ln_sla) THEN 
    253             jtype = jtype + 1 
    254             clsurffiles(jtype,:) = cn_slafbfiles(:) 
    255             cobstypessurf(jtype) = 'sla   ' 
    256             ifilessurf(jtype) = 0 
    257             DO jfile = 1, jpmaxnfiles 
    258                IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
    259                   ifilessurf(jtype) = ifilessurf(jtype) + 1 
    260             END DO 
    261          ENDIF 
    262          IF (ln_sst) THEN 
    263             jtype = jtype + 1 
    264             clsurffiles(jtype,:) = cn_sstfbfiles(:) 
    265             cobstypessurf(jtype) = 'sst   ' 
    266             ifilessurf(jtype) = 0 
    267             DO jfile = 1, jpmaxnfiles 
    268                IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
    269                   ifilessurf(jtype) = ifilessurf(jtype) + 1 
    270             END DO 
    271          ENDIF 
    272 #if defined key_lim3 
    273          IF (ln_sic) THEN 
    274             jtype = jtype + 1 
    275             clsurffiles(jtype,:) = cn_sicfbfiles(:) 
    276             cobstypessurf(jtype) = 'sic   ' 
    277             ifilessurf(jtype) = 0 
    278             DO jfile = 1, jpmaxnfiles 
    279                IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
    280                   ifilessurf(jtype) = ifilessurf(jtype) + 1 
    281             END DO 
    282          ENDIF 
    283 #endif 
    284  
    285       ENDIF 
    286  
    287       !Write namelist settings to stdout 
     233 
    288234      IF(lwp) THEN 
    289235         WRITE(numout,*) 
     
    297243         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic 
    298244         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
    299          WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ',ln_grid_global 
    300          WRITE(numout,*) '             Logical switch for SST bias correction         ln_sstbias = ', ln_sstbias  
    301          WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup 
     245         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss 
     246         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global 
     247         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 
    302248         IF (ln_grid_search_lookup) & 
    303249            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile 
     
    307253         WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint 
    308254         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
     255         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
    309256         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc 
    310257         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr 
    311258         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff 
    312259         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias 
     260         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias 
    313261         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis 
    314262         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes 
    315263         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
    316          WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
    317  
    318          IF ( nproftypes > 0 ) THEN 
    319             DO jtype = 1, nproftypes 
    320                DO jfile = 1, ifilesprof(jtype) 
    321                   WRITE(numout,'(1X,2A)') '             '//cobstypesprof(jtype)//' input observation file names  = ', & 
    322                      TRIM(clproffiles(jtype,jfile)) 
    323                END DO 
    324             END DO 
     264      ENDIF 
     265      !----------------------------------------------------------------------- 
     266      ! Set up list of observation types to be used 
     267      ! and the files associated with each type 
     268      !----------------------------------------------------------------------- 
     269 
     270      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
     271      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss /) ) 
     272 
     273      IF (ln_sstbias) THEN  
     274         lmask(:) = .FALSE.  
     275         WHERE (cn_sstbiasfiles(:) /= '') lmask(:) = .TRUE.  
     276         jnumsstbias = COUNT(lmask)  
     277         lmask(:) = .FALSE.  
     278      ENDIF       
     279 
     280      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
     281         IF(lwp) WRITE(numout,cform_war) 
     282         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
     283            &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 
     284            &                    ' are set to .FALSE. so turning off calls to dia_obs' 
     285         nwarn = nwarn + 1 
     286         ln_diaobs = .FALSE. 
     287         RETURN 
     288      ENDIF 
     289 
     290      IF ( nproftypes > 0 ) THEN 
     291 
     292         ALLOCATE( cobstypesprof(nproftypes) ) 
     293         ALLOCATE( ifilesprof(nproftypes) ) 
     294         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
     295 
     296         jtype = 0 
     297         IF (ln_t3d .OR. ln_s3d) THEN 
     298            jtype = jtype + 1 
     299            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', & 
     300               &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    325301         ENDIF 
    326  
    327          WRITE(numout,*)'          Number of surface obs types: ',nsurftypes 
    328          IF ( nsurftypes > 0 ) THEN 
    329             DO jtype = 1, nsurftypes 
    330                DO jfile = 1, ifilessurf(jtype) 
    331                   WRITE(numout,'(1X,2A)') '             '//cobstypessurf(jtype)//' input observation file names  = ', & 
    332                      TRIM(clsurffiles(jtype,jfile)) 
    333                END DO 
    334             END DO 
     302         IF (ln_vel3d) THEN 
     303            jtype = jtype + 1 
     304            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', & 
     305               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    335306         ENDIF 
    336          WRITE(numout,*) '~~~~~~~~~~~~' 
    337  
    338       ENDIF 
     307 
     308      ENDIF 
     309 
     310      IF ( nsurftypes > 0 ) THEN 
     311 
     312         ALLOCATE( cobstypessurf(nsurftypes) ) 
     313         ALLOCATE( ifilessurf(nsurftypes) ) 
     314         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
     315         ALLOCATE(n2dintsurf(nsurftypes)) 
     316         ALLOCATE(zavglamscl(nsurftypes)) 
     317         ALLOCATE(zavgphiscl(nsurftypes)) 
     318         ALLOCATE(lfpindegs(nsurftypes)) 
     319         ALLOCATE(llnightav(nsurftypes)) 
     320 
     321         jtype = 0 
     322         IF (ln_sla) THEN 
     323            jtype = jtype + 1 
     324            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', & 
     325               &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     326            CALL obs_setinterpopts( nsurftypes, jtype, 'sla   ',      & 
     327               &                  nn_2dint, nn_2dint_sla,             & 
     328               &                  rn_sla_avglamscl, rn_sla_avgphiscl, & 
     329               &                  ln_sla_fp_indegs, .FALSE.,          & 
     330               &                  n2dintsurf, zavglamscl, zavgphiscl, & 
     331               &                  lfpindegs, llnightav ) 
     332         ENDIF 
     333         IF (ln_sst) THEN 
     334            jtype = jtype + 1 
     335            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', & 
     336               &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     337            CALL obs_setinterpopts( nsurftypes, jtype, 'sst   ',      & 
     338               &                  nn_2dint, nn_2dint_sst,             & 
     339               &                  rn_sst_avglamscl, rn_sst_avgphiscl, & 
     340               &                  ln_sst_fp_indegs, ln_sstnight,      & 
     341               &                  n2dintsurf, zavglamscl, zavgphiscl, & 
     342               &                  lfpindegs, llnightav ) 
     343         ENDIF 
     344#if defined key_lim3 || defined key_cice 
     345         IF (ln_sic) THEN 
     346            jtype = jtype + 1 
     347            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', & 
     348               &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     349            CALL obs_setinterpopts( nsurftypes, jtype, 'sic   ',      & 
     350               &                  nn_2dint, nn_2dint_sic,             & 
     351               &                  rn_sic_avglamscl, rn_sic_avgphiscl, & 
     352               &                  ln_sic_fp_indegs, .FALSE.,          & 
     353               &                  n2dintsurf, zavglamscl, zavgphiscl, & 
     354               &                  lfpindegs, llnightav ) 
     355         ENDIF 
     356#endif 
     357         IF (ln_sss) THEN 
     358            jtype = jtype + 1 
     359            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', & 
     360               &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     361            CALL obs_setinterpopts( nsurftypes, jtype, 'sss   ',      & 
     362               &                  nn_2dint, nn_2dint_sss,             & 
     363               &                  rn_sss_avglamscl, rn_sss_avgphiscl, & 
     364               &                  ln_sss_fp_indegs, .FALSE.,          & 
     365               &                  n2dintsurf, zavglamscl, zavgphiscl, & 
     366               &                  lfpindegs, llnightav ) 
     367         ENDIF 
     368 
     369      ENDIF 
     370 
     371 
    339372 
    340373      !----------------------------------------------------------------------- 
     
    356389      ENDIF 
    357390 
    358       IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN 
     391      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN 
    359392         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 
    360393            &                    ' is not available') 
     
    421454               &               jpi, jpj, jpk, & 
    422455               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  & 
    423                &               ln_nea, kdailyavtypes = nn_profdavtypes ) 
     456               &               ln_nea, ln_bound_reject, & 
     457               &               kdailyavtypes = nn_profdavtypes ) 
    424458 
    425459         END DO 
     
    440474            nvarssurf(jtype) = 1 
    441475            nextrsurf(jtype) = 0 
    442             llnightav = .FALSE. 
     476            llnightav(jtype) = .FALSE. 
    443477            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 
    444             IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight 
     478            IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav(jtype) = ln_sstnight 
    445479 
    446480            !Read in surface obs types 
     
    448482               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
    449483               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
    450                &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 
    451           
    452           
    453             CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 
     484               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 
     485 
     486            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
    454487 
    455488            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
    456                CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 
    457                IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 
     489               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 
     490               IF ( ln_altbias ) & 
     491                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 
    458492            ENDIF 
     493 
    459494            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 
    460                !Read in bias field and correct SST. 
    461                IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 
    462                                                      "  but no bias"// & 
    463                                                      " files to read in")    
    464                   CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, & 
    465                                         jnumsstbias, cn_sstbias_files(1:jnumsstbias) ) 
     495               jnumsstbias = 0 
     496               DO jfile = 1, jpmaxnfiles 
     497                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 
     498                     &  jnumsstbias = jnumsstbias + 1 
     499               END DO 
     500               IF ( jnumsstbias == 0 ) THEN 
     501                  CALL ctl_stop("ln_sstbias set but no bias files to read in")     
     502               ENDIF 
     503 
     504               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &  
     505                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )  
     506 
    466507            ENDIF 
    467508         END DO 
     
    512553      USE ice    , ONLY : at_i                ! LIM3 Ice model variables 
    513554#endif 
     555#if defined key_cice 
     556      USE sbc_oce, ONLY : fr_i     ! ice fraction 
     557#endif 
     558 
    514559      IMPLICIT NONE 
    515560 
     
    528573         & zprofmask2              ! Mask associated with zprofvar2 
    529574      REAL(wp), POINTER, DIMENSION(:,:) :: & 
    530          & zsurfvar                ! Model values equivalent to surface ob. 
     575         & zsurfvar, &             ! Model values equivalent to surface ob. 
     576         & zsurfmask               ! Mask associated with surface variable 
    531577      REAL(wp), POINTER, DIMENSION(:,:) :: & 
    532578         & zglam1,    &            ! Model longitudes for prof variable 1 
     
    534580         & zgphi1,    &            ! Model latitudes for prof variable 1 
    535581         & zgphi2                  ! Model latitudes for prof variable 2 
    536 #if ! defined key_lim3 
    537       REAL(wp), POINTER, DIMENSION(:,:) :: at_i 
    538 #endif 
    539       LOGICAL :: llnightav        ! Logical for calculating night-time average 
    540582 
    541583      !Allocate local work arrays 
     
    545587      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 
    546588      CALL wrk_alloc( jpi, jpj, zsurfvar ) 
     589      CALL wrk_alloc( jpi, jpj, zsurfmask ) 
    547590      CALL wrk_alloc( jpi, jpj, zglam1 ) 
    548591      CALL wrk_alloc( jpi, jpj, zglam2 ) 
    549592      CALL wrk_alloc( jpi, jpj, zgphi1 ) 
    550593      CALL wrk_alloc( jpi, jpj, zgphi2 ) 
    551 #if ! defined key_lim3 
    552       CALL wrk_alloc(jpi,jpj,at_i)  
    553 #endif 
    554594      !----------------------------------------------------------------------- 
    555595 
     
    562602      idaystp = NINT( rday / rdt ) 
    563603 
    564       !----------------------------------------------------------------------- 
    565       ! No LIM => at_i == 0.0_wp 
    566       !----------------------------------------------------------------------- 
    567 #if ! defined key_lim3 
    568       at_i(:,:) = 0.0_wp 
    569 #endif 
    570604      !----------------------------------------------------------------------- 
    571605      ! Call the profile and surface observation operators 
     
    595629               zgphi1(:,:) = gphiu(:,:) 
    596630               zgphi2(:,:) = gphiv(:,:) 
     631            CASE DEFAULT 
     632               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 
    597633            END SELECT 
    598634 
    599             IF( ln_zco .OR. ln_zps ) THEN  
    600                CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
    601                   &               nit000, idaystp,                         & 
    602                   &               zprofvar1, zprofvar2,                    & 
    603                   &               gdept_1d, zprofmask1, zprofmask2,        & 
    604                   &               zglam1, zglam2, zgphi1, zgphi2,          & 
    605                   &               nn_1dint, nn_2dint,                      & 
    606                   &               kdailyavtypes = nn_profdavtypes ) 
    607             ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') THEN 
    608                !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION 
    609                CALL obs_pro_sco_opt( profdataqc(jtype),                    &  
    610                   &              kstp, jpi, jpj, jpk, nit000, idaystp,   &  
    611                   &              zprofvar1, zprofvar2,                   &  
    612                   &              gdept_n(:,:,:), gdepw_n(:,:,:),           & 
    613                   &              tmask, nn_1dint, nn_2dint,              &  
    614                   &              kdailyavtypes = nn_profdavtypes )  
    615             ELSE 
    616                CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// & 
    617                          'yet working for velocity data (turn off velocity observations') 
    618             ENDIF 
     635            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
     636               &               nit000, idaystp,                         & 
     637               &               zprofvar1, zprofvar2,                    & 
     638               &               gdept_n(:,:,:), gdepw_n(:,:,:),            &  
     639               &               zprofmask1, zprofmask2,                  & 
     640               &               zglam1, zglam2, zgphi1, zgphi2,          & 
     641               &               nn_1dint, nn_2dint,                      & 
     642               &               kdailyavtypes = nn_profdavtypes ) 
    619643 
    620644         END DO 
     
    625649 
    626650         DO jtype = 1, nsurftypes 
     651 
     652            !Defaults which might be changed 
     653            zsurfmask(:,:) = tmask(:,:,1) 
    627654 
    628655            SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
    629656            CASE('sst') 
    630657               zsurfvar(:,:) = tsn(:,:,1,jp_tem) 
    631                llnightav = ln_sstnight 
    632658            CASE('sla') 
    633659               zsurfvar(:,:) = sshn(:,:) 
    634                llnightav = .FALSE. 
    635 #if defined key_lim3 
     660            CASE('sss') 
     661               zsurfvar(:,:) = tsn(:,:,1,jp_sal) 
    636662            CASE('sic') 
    637663               IF ( kstp == 0 ) THEN 
     
    646672                  CYCLE 
    647673               ELSE 
    648                   zsurfvar(:,:) = at_i(:,:) 
     674#if defined key_cice 
     675                  zsurfvar(:,:) = fr_i(:,:) 
     676#elif defined key_lim2 || defined key_lim3 
     677                  zsurfvar(:,:) = 1._wp - frld(:,:) 
     678#else 
     679               CALL ctl_stop( ' Trying to run sea-ice observation operator', & 
     680                  &           ' but no sea-ice model appears to have been defined' ) 
     681#endif 
    649682               ENDIF 
    650683 
    651                llnightav = .FALSE. 
    652 #endif 
    653684            END SELECT 
    654685 
    655686            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    656                &               nit000, idaystp, zsurfvar, tmask(:,:,1), & 
    657                &               nn_2dint, llnightav ) 
     687               &               nit000, idaystp, zsurfvar, zsurfmask,    & 
     688               &               n2dintsurf(jtype), llnightav(jtype),     & 
     689               &               zavglamscl(jtype), zavgphiscl(jtype),     & 
     690               &               lfpindegs(jtype) ) 
    658691 
    659692         END DO 
     
    666699      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 
    667700      CALL wrk_dealloc( jpi, jpj, zsurfvar ) 
     701      CALL wrk_dealloc( jpi, jpj, zsurfmask ) 
    668702      CALL wrk_dealloc( jpi, jpj, zglam1 ) 
    669703      CALL wrk_dealloc( jpi, jpj, zglam2 ) 
    670704      CALL wrk_dealloc( jpi, jpj, zgphi1 ) 
    671705      CALL wrk_dealloc( jpi, jpj, zgphi2 ) 
    672 #if ! defined key_lim3 
    673       CALL wrk_dealloc(jpi,jpj,at_i) 
    674 #endif 
    675706 
    676707   END SUBROUTINE dia_obs 
     
    789820 
    790821      IF ( nsurftypes > 0 ) & 
    791          &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 
     822         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 
     823         &               n2dintsurf, zavglamscl, zavgphiscl, lfpindegs, llnightav ) 
    792824 
    793825   END SUBROUTINE dia_obs_dealloc 
     
    938970   END SUBROUTINE fin_date 
    939971    
     972    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, & 
     973       &                         cfilestype, ifiles, cobstypes, cfiles ) 
     974 
     975    INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types 
     976    INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 
     977    INTEGER, INTENT(IN) :: jtype       ! Index of the current type of obs 
     978    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
     979       &                   ifiles      ! Out appended number of files for this type 
     980 
     981    CHARACTER(len=6), INTENT(IN) :: ctypein  
     982    CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: & 
     983       &                   cfilestype  ! In list of files for this obs type 
     984    CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: & 
     985       &                   cobstypes   ! Out appended list of obs types 
     986    CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: & 
     987       &                   cfiles      ! Out appended list of files for all types 
     988 
     989    !Local variables 
     990    INTEGER :: jfile 
     991 
     992    cfiles(jtype,:) = cfilestype(:) 
     993    cobstypes(jtype) = ctypein 
     994    ifiles(jtype) = 0 
     995    DO jfile = 1, jpmaxnfiles 
     996       IF ( trim(cfiles(jtype,jfile)) /= '' ) & 
     997                 ifiles(jtype) = ifiles(jtype) + 1 
     998    END DO 
     999 
     1000    IF ( ifiles(jtype) == 0 ) THEN 
     1001         CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)//   & 
     1002            &           ' set to true but no files available to read' ) 
     1003    ENDIF 
     1004 
     1005    IF(lwp) THEN     
     1006       WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:' 
     1007       DO jfile = 1, ifiles(jtype) 
     1008          WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile)) 
     1009       END DO 
     1010    ENDIF 
     1011 
     1012    END SUBROUTINE obs_settypefiles 
     1013 
     1014    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             & 
     1015               &                  n2dint_default, n2dint_type,        & 
     1016               &                  zavglamscl_type, zavgphiscl_type,   & 
     1017               &                  lfp_indegs_type, lavnight_type,     & 
     1018               &                  n2dint, zavglamscl, zavgphiscl,     & 
     1019               &                  lfpindegs, lavnight ) 
     1020 
     1021    INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types 
     1022    INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs 
     1023    INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type 
     1024    INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type 
     1025    REAL(wp), INTENT(IN) :: & 
     1026       &                    zavglamscl_type, & !E/W diameter of obs footprint for this type 
     1027       &                    zavgphiscl_type    !N/S diameter of obs footprint for this type 
     1028    LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres 
     1029    LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average 
     1030    CHARACTER(len=6), INTENT(IN) :: ctypein  
     1031 
     1032    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1033       &                    n2dint  
     1034    REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 
     1035       &                    zavglamscl, zavgphiscl 
     1036    LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1037       &                    lfpindegs, lavnight 
     1038 
     1039    lavnight(jtype) = lavnight_type 
     1040 
     1041    IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN 
     1042       n2dint(jtype) = n2dint_type 
     1043    ELSE 
     1044       n2dint(jtype) = n2dint_default 
     1045    ENDIF 
     1046 
     1047    ! For averaging observation footprints set options for size of footprint  
     1048    IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 
     1049       IF ( zavglamscl_type > 0._wp ) THEN 
     1050          zavglamscl(jtype) = zavglamscl_type 
     1051       ELSE 
     1052          CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1053                         'scale (zavglamscl) for observation type '//TRIM(ctypein) )       
     1054       ENDIF 
     1055 
     1056       IF ( zavgphiscl_type > 0._wp ) THEN 
     1057          zavgphiscl(jtype) = zavgphiscl_type 
     1058       ELSE 
     1059          CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1060                         'scale (zavgphiscl) for observation type '//TRIM(ctypein) )       
     1061       ENDIF 
     1062 
     1063       lfpindegs(jtype) = lfp_indegs_type  
     1064 
     1065    ENDIF 
     1066 
     1067    ! Write out info  
     1068    IF(lwp) THEN 
     1069       IF ( n2dint(jtype) <= 4 ) THEN 
     1070          WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1071             &            ' model counterparts will be interpolated horizontally' 
     1072       ELSE IF ( n2dint(jtype) <= 6 ) THEN 
     1073          WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1074             &            ' model counterparts will be averaged horizontally' 
     1075          WRITE(numout,*) '             '//'    with E/W scale: ',zavglamscl(jtype) 
     1076          WRITE(numout,*) '             '//'    with N/S scale: ',zavgphiscl(jtype) 
     1077          IF ( lfpindegs(jtype) ) THEN 
     1078              WRITE(numout,*) '             '//'    (in degrees)' 
     1079          ELSE 
     1080              WRITE(numout,*) '             '//'    (in metres)' 
     1081          ENDIF 
     1082       ENDIF 
     1083    ENDIF 
     1084 
     1085    END SUBROUTINE obs_setinterpopts 
     1086 
    9401087END MODULE diaobs 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_mpp.F90

    r6140 r9023  
    142142      !! 
    143143 
    144       iobsp=kobsp 
     144      iobsp(:)=kobsp(:) 
    145145 
    146146      WHERE( iobsp(:) == -1 ) 
     
    148148      END WHERE 
    149149 
    150       iobsp=-1*iobsp 
     150      iobsp(:)=-1*iobsp(:) 
    151151 
    152152      CALL obs_mpp_max_integer( iobsp, kno ) 
    153153 
    154       kobsp=-1*iobsp 
     154      kobsp(:)=-1*iobsp(:) 
    155155 
    156156      isum=0 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r7646 r9023  
    99   !!   obs_prof_opt :    Compute the model counterpart of profile data 
    1010   !!   obs_surf_opt :    Compute the model counterpart of surface data 
    11    !!   obs_pro_sco_opt: Compute the model counterpart of temperature and  
    12    !!                    salinity observations from profiles in generalised  
    13    !!                    vertical coordinates  
    1411   !!---------------------------------------------------------------------- 
    1512 
     
    2219      & obs_int_h2d, & 
    2320      & obs_int_h2d_init 
     21   USE obs_averg_h2d, ONLY : &    ! Horizontal averaging to the obs footprint 
     22      & obs_avg_h2d, & 
     23      & obs_avg_h2d_init, & 
     24      & obs_max_fpsize 
    2425   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the obs pt 
    2526      & obs_int_z1d,    & 
    2627      & obs_int_z1d_spl 
    27    USE obs_const,  ONLY :     & 
    28       & obfillflt      ! Fillvalue    
     28   USE obs_const,  ONLY :    &    ! Obs fill value 
     29      & obfillflt 
    2930   USE dom_oce,       ONLY : & 
    30       & glamt, glamu, glamv, & 
    31       & gphit, gphiu, gphiv, &  
    32       & gdept_n, gdept_0  
    33    USE lib_mpp,       ONLY : & 
     31      & glamt, glamf, & 
     32      & gphit, gphif 
     33   USE lib_mpp,       ONLY : &    ! Warning and stopping routines 
    3434      & ctl_warn, ctl_stop 
     35   USE sbcdcy,        ONLY : &    ! For calculation of where it is night-time 
     36      & sbc_dcy, nday_qsr 
    3537   USE obs_grid,      ONLY : &  
    3638      & obs_level_search      
    37    USE sbcdcy,        ONLY : &    ! For calculation of where it is night-time 
    38       & sbc_dcy, nday_qsr 
    3939 
    4040   IMPLICIT NONE 
     
    4444 
    4545   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs 
    46       &   obs_pro_sco_opt, &  ! Compute the model counterpart of profile observations  
    4746      &   obs_surf_opt     ! Compute the model counterpart of surface obs 
    4847 
     
    5857CONTAINS 
    5958 
     59 
    6060   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          & 
    6161      &                     kit000, kdaystp,                      & 
    62       &                     pvar1, pvar2, pgdept, pmask1, pmask2, & 
     62      &                     pvar1, pvar2, pgdept, pgdepw,         & 
     63      &                     pmask1, pmask2,                       &   
    6364      &                     plam1, plam2, pphi1, pphi2,           & 
    6465      &                     k1dint, k2dint, kdailyavtypes ) 
     
    111112      !!      ! 07-03 (K. Mogensen) General handling of profiles 
    112113      !!      ! 15-02 (M. Martin) Combined routine for all profile types 
     114      !!      ! 17-02 (M. Martin) Include generalised vertical coordinate changes 
    113115      !!----------------------------------------------------------------------- 
    114116 
     
    140142         & pphi1,    &               ! Model latitudes for variable 1 
    141143         & pphi2                     ! Model latitudes for variable 2 
    142       REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 
    143          & pgdept                    ! Model array of depth levels 
     144      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
     145         & pgdept, &                 ! Model array of depth T levels  
     146         & pgdepw                    ! Model array of depth W levels  
    144147      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    145148         & kdailyavtypes             ! Types for daily averages 
     
    156159      INTEGER ::   iend 
    157160      INTEGER ::   iobs 
     161      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes  
     162      INTEGER ::   inum_obs 
    158163      INTEGER, DIMENSION(imaxavtypes) :: & 
    159164         & idailyavtypes 
     
    163168         & igrdj1, & 
    164169         & igrdj2 
     170      INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 
     171 
    165172      REAL(KIND=wp) :: zlam 
    166173      REAL(KIND=wp) :: zphi 
     
    171178         & zobsk,    & 
    172179         & zobs2k 
    173       REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 
     180      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
    174181         & zweig1, & 
    175          & zweig2 
     182         & zweig2, & 
     183         & zweig 
    176184      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    177185         & zmask1, & 
    178186         & zmask2, & 
    179          & zint1, & 
    180          & zint2, & 
    181          & zinm1, & 
    182          & zinm2 
     187         & zint1,  & 
     188         & zint2,  & 
     189         & zinm1,  & 
     190         & zinm2,  & 
     191         & zgdept, &  
     192         & zgdepw 
    183193      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    184194         & zglam1, & 
     
    186196         & zgphi1, & 
    187197         & zgphi2 
     198      REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2    
     199      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 
     200 
    188201      LOGICAL :: ld_dailyav 
    189202 
     
    266279         & zmask1(2,2,kpk,ipro),  & 
    267280         & zmask2(2,2,kpk,ipro),  & 
    268          & zint1(2,2,kpk,ipro),  & 
    269          & zint2(2,2,kpk,ipro)   & 
     281         & zint1(2,2,kpk,ipro),   & 
     282         & zint2(2,2,kpk,ipro),   & 
     283         & zgdept(2,2,kpk,ipro),  &  
     284         & zgdepw(2,2,kpk,ipro)   &  
    270285         & ) 
    271286 
     
    290305      END DO 
    291306 
     307      ! Initialise depth arrays 
     308      zgdept(:,:,:,:) = 0.0 
     309      zgdepw(:,:,:,:) = 0.0 
     310 
    292311      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 
    293312      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 
     
    300319      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
    301320 
     321      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept )  
     322      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw )  
     323 
    302324      ! At the end of the day also get interpolated means 
    303325      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
     
    314336 
    315337      ENDIF 
     338 
     339      ! Return if no observations to process  
     340      ! Has to be done after comm commands to ensure processors  
     341      ! stay in sync  
     342      IF ( ipro == 0 ) RETURN  
    316343 
    317344      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
     
    339366         zphi = prodatqc%rphi(jobs) 
    340367 
    341          ! Horizontal weights and vertical mask 
    342  
     368         ! Horizontal weights  
     369         ! Masked values are calculated later.   
    343370         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    344371 
    345             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     372            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    346373               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), & 
    347                &                   zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 
     374               &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 
    348375 
    349376         ENDIF 
     
    351378         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    352379 
    353             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     380            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    354381               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
    355                &                   zmask2(:,:,:,iobs), zweig2, zobsmask2 ) 
     382               &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 
    356383  
    357384         ENDIF 
     
    365392               IF ( idayend == 0 )  THEN 
    366393                  ! Daily averaged data 
    367                   CALL obs_int_h2d( kpk, kpk,      & 
    368                      &              zweig1, zinm1(:,:,:,iobs), zobsk ) 
    369  
    370                ENDIF 
    371  
    372             ELSE  
    373  
    374                ! Point data 
    375                CALL obs_int_h2d( kpk, kpk,      & 
    376                   &              zweig1, zint1(:,:,:,iobs), zobsk ) 
    377  
    378             ENDIF 
    379  
    380             !------------------------------------------------------------- 
    381             ! Compute vertical second-derivative of the interpolating  
    382             ! polynomial at obs points 
    383             !------------------------------------------------------------- 
    384  
    385             IF ( k1dint == 1 ) THEN 
    386                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   & 
    387                   &                  pgdept, zobsmask1 ) 
    388             ENDIF 
    389  
    390             !----------------------------------------------------------------- 
    391             !  Vertical interpolation to the observation point 
    392             !----------------------------------------------------------------- 
    393             ista = prodatqc%npvsta(jobs,1) 
    394             iend = prodatqc%npvend(jobs,1) 
    395             CALL obs_int_z1d( kpk,                & 
    396                & prodatqc%var(1)%mvk(ista:iend),  & 
    397                & k1dint, iend - ista + 1,         & 
    398                & prodatqc%var(1)%vdep(ista:iend), & 
    399                & zobsk, zobs2k,                   & 
    400                & prodatqc%var(1)%vmod(ista:iend), & 
    401                & pgdept, zobsmask1 ) 
    402  
    403          ENDIF 
    404  
     394 
     395                  ! vertically interpolate all 4 corners  
     396                  ista = prodatqc%npvsta(jobs,1)  
     397                  iend = prodatqc%npvend(jobs,1)  
     398                  inum_obs = iend - ista + 1  
     399                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     400 
     401                  DO iin=1,2  
     402                     DO ijn=1,2  
     403 
     404                        IF ( k1dint == 1 ) THEN  
     405                           CALL obs_int_z1d_spl( kpk, &  
     406                              &     zinm1(iin,ijn,:,iobs), &  
     407                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     408                              &     zmask1(iin,ijn,:,iobs))  
     409                        ENDIF  
     410        
     411                        CALL obs_level_search(kpk, &  
     412                           &    zgdept(iin,ijn,:,iobs), &  
     413                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     414                           &    iv_indic)  
     415 
     416                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     417                           &    prodatqc%var(1)%vdep(ista:iend), &  
     418                           &    zinm1(iin,ijn,:,iobs), &  
     419                           &    zobs2k, interp_corner(iin,ijn,:), &  
     420                           &    zgdept(iin,ijn,:,iobs), &  
     421                           &    zmask1(iin,ijn,:,iobs))  
     422        
     423                     ENDDO  
     424                  ENDDO  
     425 
     426               ENDIF !idayend 
     427 
     428            ELSE    
     429 
     430               ! Point data  
     431      
     432               ! vertically interpolate all 4 corners  
     433               ista = prodatqc%npvsta(jobs,1)  
     434               iend = prodatqc%npvend(jobs,1)  
     435               inum_obs = iend - ista + 1  
     436               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     437               DO iin=1,2   
     438                  DO ijn=1,2  
     439                     
     440                     IF ( k1dint == 1 ) THEN  
     441                        CALL obs_int_z1d_spl( kpk, &  
     442                           &    zint1(iin,ijn,:,iobs),&  
     443                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     444                           &    zmask1(iin,ijn,:,iobs))  
     445   
     446                     ENDIF  
     447        
     448                     CALL obs_level_search(kpk, &  
     449                         &        zgdept(iin,ijn,:,iobs),&  
     450                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     451                         &        iv_indic)  
     452 
     453                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     454                         &          prodatqc%var(1)%vdep(ista:iend),     &  
     455                         &          zint1(iin,ijn,:,iobs),            &  
     456                         &          zobs2k,interp_corner(iin,ijn,:), &  
     457                         &          zgdept(iin,ijn,:,iobs),         &  
     458                         &          zmask1(iin,ijn,:,iobs) )       
     459          
     460                  ENDDO  
     461               ENDDO  
     462              
     463            ENDIF  
     464 
     465            !-------------------------------------------------------------  
     466            ! Compute the horizontal interpolation for every profile level  
     467            !-------------------------------------------------------------  
     468              
     469            DO ikn=1,inum_obs  
     470               iend=ista+ikn-1 
     471                   
     472               zweig(:,:,1) = 0._wp  
     473    
     474               ! This code forces the horizontal weights to be   
     475               ! zero IF the observation is below the bottom of the   
     476               ! corners of the interpolation nodes, Or if it is in   
     477               ! the mask. This is important for observations near   
     478               ! steep bathymetry  
     479               DO iin=1,2  
     480                  DO ijn=1,2  
     481      
     482                     depth_loop1: DO ik=kpk,2,-1  
     483                        IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     484                             
     485                           zweig(iin,ijn,1) = &   
     486                              & zweig1(iin,ijn,1) * &  
     487                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     488                              &  - prodatqc%var(1)%vdep(iend)),0._wp)  
     489                             
     490                           EXIT depth_loop1  
     491 
     492                        ENDIF  
     493 
     494                     ENDDO depth_loop1  
     495      
     496                  ENDDO  
     497               ENDDO  
     498    
     499               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     500                  &              prodatqc%var(1)%vmod(iend:iend) )  
     501 
     502                  ! Set QC flag for any observations found below the bottom 
     503                  ! needed as the check here is more strict than that in obs_prep 
     504               IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4 
     505  
     506            ENDDO  
     507  
     508            DEALLOCATE(interp_corner,iv_indic)  
     509           
     510         ENDIF  
     511 
     512         ! For the second variable 
    405513         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    406514 
     
    410518 
    411519               IF ( idayend == 0 )  THEN 
    412  
    413520                  ! Daily averaged data 
    414                   CALL obs_int_h2d( kpk, kpk,      & 
    415                      &              zweig2, zinm2(:,:,:,iobs), zobsk ) 
    416  
    417                ENDIF 
    418  
    419             ELSE 
    420  
    421                ! Point data 
    422                CALL obs_int_h2d( kpk, kpk,      & 
    423                   &              zweig2, zint2(:,:,:,iobs), zobsk ) 
    424  
    425             ENDIF 
    426  
    427  
    428             !------------------------------------------------------------- 
    429             ! Compute vertical second-derivative of the interpolating  
    430             ! polynomial at obs points 
    431             !------------------------------------------------------------- 
    432  
    433             IF ( k1dint == 1 ) THEN 
    434                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 
    435                   &                  pgdept, zobsmask2 ) 
    436             ENDIF 
    437  
    438             !---------------------------------------------------------------- 
    439             !  Vertical interpolation to the observation point 
    440             !---------------------------------------------------------------- 
    441             ista = prodatqc%npvsta(jobs,2) 
    442             iend = prodatqc%npvend(jobs,2) 
    443             CALL obs_int_z1d( kpk, & 
    444                & prodatqc%var(2)%mvk(ista:iend),& 
    445                & k1dint, iend - ista + 1, & 
    446                & prodatqc%var(2)%vdep(ista:iend),& 
    447                & zobsk, zobs2k, & 
    448                & prodatqc%var(2)%vmod(ista:iend),& 
    449                & pgdept, zobsmask2 ) 
    450  
    451          ENDIF 
    452  
    453       END DO 
     521 
     522                  ! vertically interpolate all 4 corners  
     523                  ista = prodatqc%npvsta(jobs,2)  
     524                  iend = prodatqc%npvend(jobs,2)  
     525                  inum_obs = iend - ista + 1  
     526                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     527 
     528                  DO iin=1,2  
     529                     DO ijn=1,2  
     530 
     531                        IF ( k1dint == 1 ) THEN  
     532                           CALL obs_int_z1d_spl( kpk, &  
     533                              &     zinm2(iin,ijn,:,iobs), &  
     534                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     535                              &     zmask2(iin,ijn,:,iobs))  
     536                        ENDIF  
     537        
     538                        CALL obs_level_search(kpk, &  
     539                           &    zgdept(iin,ijn,:,iobs), &  
     540                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     541                           &    iv_indic)  
     542 
     543                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     544                           &    prodatqc%var(2)%vdep(ista:iend), &  
     545                           &    zinm2(iin,ijn,:,iobs), &  
     546                           &    zobs2k, interp_corner(iin,ijn,:), &  
     547                           &    zgdept(iin,ijn,:,iobs), &  
     548                           &    zmask2(iin,ijn,:,iobs))  
     549        
     550                     ENDDO  
     551                  ENDDO  
     552 
     553               ENDIF !idayend 
     554 
     555            ELSE    
     556 
     557               ! Point data  
     558      
     559               ! vertically interpolate all 4 corners  
     560               ista = prodatqc%npvsta(jobs,2)  
     561               iend = prodatqc%npvend(jobs,2)  
     562               inum_obs = iend - ista + 1  
     563               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     564               DO iin=1,2   
     565                  DO ijn=1,2  
     566                     
     567                     IF ( k1dint == 1 ) THEN  
     568                        CALL obs_int_z1d_spl( kpk, &  
     569                           &    zint2(iin,ijn,:,iobs),&  
     570                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     571                           &    zmask2(iin,ijn,:,iobs))  
     572   
     573                     ENDIF  
     574        
     575                     CALL obs_level_search(kpk, &  
     576                         &        zgdept(iin,ijn,:,iobs),&  
     577                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     578                         &        iv_indic)  
     579 
     580                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     581                         &          prodatqc%var(2)%vdep(ista:iend),     &  
     582                         &          zint2(iin,ijn,:,iobs),            &  
     583                         &          zobs2k,interp_corner(iin,ijn,:), &  
     584                         &          zgdept(iin,ijn,:,iobs),         &  
     585                         &          zmask2(iin,ijn,:,iobs) )       
     586          
     587                  ENDDO  
     588               ENDDO  
     589              
     590            ENDIF  
     591 
     592            !-------------------------------------------------------------  
     593            ! Compute the horizontal interpolation for every profile level  
     594            !-------------------------------------------------------------  
     595              
     596            DO ikn=1,inum_obs  
     597               iend=ista+ikn-1 
     598                   
     599               zweig(:,:,1) = 0._wp  
     600    
     601               ! This code forces the horizontal weights to be   
     602               ! zero IF the observation is below the bottom of the   
     603               ! corners of the interpolation nodes, Or if it is in   
     604               ! the mask. This is important for observations near   
     605               ! steep bathymetry  
     606               DO iin=1,2  
     607                  DO ijn=1,2  
     608      
     609                     depth_loop2: DO ik=kpk,2,-1  
     610                        IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     611                             
     612                           zweig(iin,ijn,1) = &   
     613                              & zweig2(iin,ijn,1) * &  
     614                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     615                              &  - prodatqc%var(2)%vdep(iend)),0._wp)  
     616                             
     617                           EXIT depth_loop2  
     618 
     619                        ENDIF  
     620 
     621                     ENDDO depth_loop2  
     622      
     623                  ENDDO  
     624               ENDDO  
     625    
     626               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     627                  &              prodatqc%var(2)%vmod(iend:iend) )  
     628 
     629                  ! Set QC flag for any observations found below the bottom 
     630                  ! needed as the check here is more strict than that in obs_prep 
     631               IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 
     632  
     633            ENDDO  
     634  
     635            DEALLOCATE(interp_corner,iv_indic)  
     636           
     637         ENDIF  
     638 
     639      ENDDO 
    454640 
    455641      ! Deallocate the data for interpolation 
     
    466652         & zmask2, & 
    467653         & zint1,  & 
    468          & zint2   & 
     654         & zint2,  & 
     655         & zgdept, & 
     656         & zgdepw  & 
    469657         & ) 
    470658 
     
    481669   END SUBROUTINE obs_prof_opt 
    482670 
    483    SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &  
    484       &                    ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, &  
    485       &                    kdailyavtypes )  
    486       !!-----------------------------------------------------------------------  
    487       !!  
    488       !!                     ***  ROUTINE obs_pro_opt  ***  
    489       !!  
    490       !! ** Purpose : Compute the model counterpart of profiles  
    491       !!              data by interpolating from the model grid to the   
    492       !!              observation point. Generalised vertical coordinate version  
    493       !!  
    494       !! ** Method  : Linearly interpolate to each observation point using   
    495       !!              the model values at the corners of the surrounding grid box.  
    496       !!  
    497       !!          First, model values on the model grid are interpolated vertically to the  
    498       !!          Depths of the profile observations.  Two vertical interpolation schemes are  
    499       !!          available:  
    500       !!          - linear       (k1dint = 0)  
    501       !!          - Cubic spline (k1dint = 1)     
    502       !!  
    503       !!  
    504       !!         Secondly the interpolated values are interpolated horizontally to the   
    505       !!         obs (lon, lat) point.  
    506       !!         Several horizontal interpolation schemes are available:  
    507       !!        - distance-weighted (great circle) (k2dint = 0)  
    508       !!        - distance-weighted (small angle)  (k2dint = 1)  
    509       !!        - bilinear (geographical grid)     (k2dint = 2)  
    510       !!        - bilinear (quadrilateral grid)    (k2dint = 3)  
    511       !!        - polynomial (quadrilateral grid)  (k2dint = 4)  
    512       !!  
    513       !!    For the cubic spline the 2nd derivative of the interpolating   
    514       !!    polynomial is computed before entering the vertical interpolation   
    515       !!    routine.  
    516       !!  
    517       !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is  
    518       !!    a daily mean model temperature field. So, we first compute  
    519       !!    the mean, then interpolate only at the end of the day.  
    520       !!  
    521       !!    This is the procedure to be used with generalised vertical model   
    522       !!    coordinates (ie s-coordinates. It is ~4x slower than the equivalent  
    523       !!    horizontal then vertical interpolation algorithm, but can deal with situations  
    524       !!    where the model levels are not flat.  
    525       !!    ONLY PERFORMED if ln_sco=.TRUE.   
    526       !!        
    527       !!    Note: the in situ temperature observations must be converted  
    528       !!    to potential temperature (the model variable) prior to  
    529       !!    assimilation.   
    530       !!??????????????????????????????????????????????????????????????  
    531       !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???  
    532       !!??????????????????????????????????????????????????????????????  
    533       !!  
    534       !! ** Action  :  
    535       !!  
    536       !! History :  
    537       !!      ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised  
    538       !!                           vertical coordinates 
    539       !!-----------------------------------------------------------------------  
    540     
    541       !! * Modules used  
    542       USE obs_profiles_def   ! Definition of storage space for profile obs.  
    543         
    544       IMPLICIT NONE  
    545   
    546       !! * Arguments  
    547       TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening  
    548       INTEGER, INTENT(IN) :: kt        ! Time step  
    549       INTEGER, INTENT(IN) :: kpi       ! Model grid parameters  
    550       INTEGER, INTENT(IN) :: kpj  
    551       INTEGER, INTENT(IN) :: kpk  
    552       INTEGER, INTENT(IN) :: kit000    ! Number of the first time step   
    553                                        !   (kit000-1 = restart time)  
    554       INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)  
    555       INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)  
    556       INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                      
    557       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
    558          & ptn,    &    ! Model temperature field  
    559          & psn,    &    ! Model salinity field  
    560          & ptmask       ! Land-sea mask  
    561       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
    562          & pgdept,  &    ! Model array of depth T levels     
    563          & pgdepw       ! Model array of depth W levels  
    564       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &  
    565          & kdailyavtypes   ! Types for daily averages  
    566        
    567       !! * Local declarations  
    568       INTEGER ::   ji  
    569       INTEGER ::   jj  
    570       INTEGER ::   jk  
    571       INTEGER ::   iico, ijco  
    572       INTEGER ::   jobs  
    573       INTEGER ::   inrc  
    574       INTEGER ::   ipro  
    575       INTEGER ::   idayend  
    576       INTEGER ::   ista  
    577       INTEGER ::   iend  
    578       INTEGER ::   iobs  
    579       INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes  
    580       INTEGER, DIMENSION(imaxavtypes) :: &  
    581          & idailyavtypes  
    582       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &  
    583          & igrdi, &  
    584          & igrdj  
    585       INTEGER :: &  
    586          & inum_obs 
    587       INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic     
    588       REAL(KIND=wp) :: zlam  
    589       REAL(KIND=wp) :: zphi  
    590       REAL(KIND=wp) :: zdaystp  
    591       REAL(KIND=wp), DIMENSION(kpk) :: &  
    592          & zobsmask, &  
    593          & zobsk,    &  
    594          & zobs2k  
    595       REAL(KIND=wp), DIMENSION(2,2,1) :: &  
    596          & zweig, &  
    597          & l_zweig  
    598       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &  
    599          & zmask, &  
    600          & zintt, &  
    601          & zints, &  
    602          & zinmt, &  
    603          & zgdept,&  
    604          & zgdepw,&  
    605          & zinms  
    606       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &  
    607          & zglam, &  
    608          & zgphi     
    609       REAL(KIND=wp), DIMENSION(1) :: zmsk_1        
    610       REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner        
    611   
    612       !------------------------------------------------------------------------  
    613       ! Local initialization   
    614       !------------------------------------------------------------------------  
    615       ! ... Record and data counters  
    616       inrc = kt - kit000 + 2  
    617       ipro = prodatqc%npstp(inrc)  
    618    
    619       ! Daily average types  
    620       IF ( PRESENT(kdailyavtypes) ) THEN  
    621          idailyavtypes(:) = kdailyavtypes(:)  
    622       ELSE  
    623          idailyavtypes(:) = -1  
    624       ENDIF  
    625   
    626       ! Initialize daily mean for first time-step  
    627       idayend = MOD( kt - kit000 + 1, kdaystp )  
    628   
    629       ! Added kt == 0 test to catch restart case   
    630       IF ( idayend == 1 .OR. kt == 0) THEN  
    631            
    632          IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt  
    633          DO jk = 1, jpk  
    634             DO jj = 1, jpj  
    635                DO ji = 1, jpi  
    636                   prodatqc%vdmean(ji,jj,jk,1) = 0.0  
    637                   prodatqc%vdmean(ji,jj,jk,2) = 0.0  
    638                END DO  
    639             END DO  
    640          END DO  
    641         
    642       ENDIF  
    643         
    644       DO jk = 1, jpk  
    645          DO jj = 1, jpj  
    646             DO ji = 1, jpi  
    647                ! Increment the temperature field for computing daily mean  
    648                prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &  
    649                &                        + ptn(ji,jj,jk)  
    650                ! Increment the salinity field for computing daily mean  
    651                prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &  
    652                &                        + psn(ji,jj,jk)  
    653             END DO  
    654          END DO  
    655       END DO  
    656      
    657       ! Compute the daily mean at the end of day  
    658       zdaystp = 1.0 / REAL( kdaystp )  
    659       IF ( idayend == 0 ) THEN  
    660          DO jk = 1, jpk  
    661             DO jj = 1, jpj  
    662                DO ji = 1, jpi  
    663                   prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &  
    664                   &                        * zdaystp  
    665                   prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &  
    666                   &                           * zdaystp  
    667                END DO  
    668             END DO  
    669          END DO  
    670       ENDIF  
    671   
    672       ! Get the data for interpolation  
    673       ALLOCATE( &  
    674          & igrdi(2,2,ipro),      &  
    675          & igrdj(2,2,ipro),      &  
    676          & zglam(2,2,ipro),      &  
    677          & zgphi(2,2,ipro),      &  
    678          & zmask(2,2,kpk,ipro),  &  
    679          & zintt(2,2,kpk,ipro),  &  
    680          & zints(2,2,kpk,ipro),  &  
    681          & zgdept(2,2,kpk,ipro), &  
    682          & zgdepw(2,2,kpk,ipro)  &  
    683          & )  
    684   
    685       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro  
    686          iobs = jobs - prodatqc%nprofup  
    687          igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1  
    688          igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1  
    689          igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1  
    690          igrdj(1,2,iobs) = prodatqc%mj(jobs,1)  
    691          igrdi(2,1,iobs) = prodatqc%mi(jobs,1)  
    692          igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1  
    693          igrdi(2,2,iobs) = prodatqc%mi(jobs,1)  
    694          igrdj(2,2,iobs) = prodatqc%mj(jobs,1)  
    695       END DO  
    696       
    697       ! Initialise depth arrays 
    698       zgdept = 0.0 
    699       zgdepw = 0.0 
    700   
    701       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, glamt, zglam )  
    702       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, gphit, zgphi )  
    703       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptmask,zmask )  
    704       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptn,   zintt )  
    705       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, psn,   zints )  
    706       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept(:,:,:), &  
    707         &                     zgdept )  
    708       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw(:,:,:), &  
    709         &                     zgdepw )  
    710   
    711       ! At the end of the day also get interpolated means  
    712       IF ( idayend == 0 ) THEN  
    713   
    714          ALLOCATE( &  
    715             & zinmt(2,2,kpk,ipro),  &  
    716             & zinms(2,2,kpk,ipro)   &  
    717             & )  
    718   
    719          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &  
    720             &                  prodatqc%vdmean(:,:,:,1), zinmt )  
    721          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &  
    722             &                  prodatqc%vdmean(:,:,:,2), zinms )  
    723   
    724       ENDIF  
    725         
    726       ! Return if no observations to process  
    727       ! Has to be done after comm commands to ensure processors  
    728       ! stay in sync  
    729       IF ( ipro == 0 ) RETURN  
    730   
    731       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro  
    732      
    733          iobs = jobs - prodatqc%nprofup  
    734      
    735          IF ( kt /= prodatqc%mstp(jobs) ) THEN  
    736               
    737             IF(lwp) THEN  
    738                WRITE(numout,*)  
    739                WRITE(numout,*) ' E R R O R : Observation',              &  
    740                   &            ' time step is not consistent with the', &  
    741                   &            ' model time step'  
    742                WRITE(numout,*) ' ========='  
    743                WRITE(numout,*)  
    744                WRITE(numout,*) ' Record  = ', jobs,                    &  
    745                   &            ' kt      = ', kt,                      &  
    746                   &            ' mstp    = ', prodatqc%mstp(jobs), &  
    747                   &            ' ntyp    = ', prodatqc%ntyp(jobs)  
    748             ENDIF  
    749             CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )  
    750          ENDIF  
    751            
    752          zlam = prodatqc%rlam(jobs)  
    753          zphi = prodatqc%rphi(jobs)  
    754            
    755          ! Horizontal weights  
    756          ! Only calculated once, for both T and S.  
    757          ! Masked values are calculated later.   
    758   
    759          IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &  
    760             & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN  
    761   
    762             CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &  
    763                &                   zglam(:,:,iobs), zgphi(:,:,iobs), &  
    764                &                   zmask(:,:,1,iobs), zweig, zmsk_1 )  
    765   
    766          ENDIF  
    767           
    768          ! IF zmsk_1 = 0; then ob is on land  
    769          IF (zmsk_1(1) < 0.1) THEN  
    770             WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask'  
    771     
    772          ELSE   
    773               
    774             ! Temperature  
    775               
    776             IF ( prodatqc%npvend(jobs,1) > 0 ) THEN   
    777      
    778                zobsk(:) = obfillflt  
    779      
    780                IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN  
    781      
    782                   IF ( idayend == 0 )  THEN  
    783                     
    784                      ! Daily averaged moored buoy (MRB) data  
    785                     
    786                      ! vertically interpolate all 4 corners  
    787                      ista = prodatqc%npvsta(jobs,1)  
    788                      iend = prodatqc%npvend(jobs,1)  
    789                      inum_obs = iend - ista + 1  
    790                      ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
    791        
    792                      DO iin=1,2  
    793                         DO ijn=1,2  
    794                                         
    795                                         
    796             
    797                            IF ( k1dint == 1 ) THEN  
    798                               CALL obs_int_z1d_spl( kpk, &  
    799                                  &     zinmt(iin,ijn,:,iobs), &  
    800                                  &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    801                                  &     zmask(iin,ijn,:,iobs))  
    802                            ENDIF  
    803         
    804                            CALL obs_level_search(kpk, &  
    805                               &    zgdept(iin,ijn,:,iobs), &  
    806                               &    inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
    807                               &    iv_indic)  
    808                            CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
    809                               &    prodatqc%var(1)%vdep(ista:iend), &  
    810                               &    zinmt(iin,ijn,:,iobs), &  
    811                               &    zobs2k, interp_corner(iin,ijn,:), &  
    812                               &    zgdept(iin,ijn,:,iobs), &  
    813                               &    zmask(iin,ijn,:,iobs))  
    814         
    815                         ENDDO  
    816                      ENDDO  
    817                     
    818                     
    819                   ELSE  
    820                  
    821                      CALL ctl_stop( ' A nonzero' //     &  
    822                         &           ' number of profile T BUOY data should' // &  
    823                         &           ' only occur at the end of a given day' )  
    824      
    825                   ENDIF  
    826           
    827                ELSE   
    828                  
    829                   ! Point data  
    830       
    831                   ! vertically interpolate all 4 corners  
    832                   ista = prodatqc%npvsta(jobs,1)  
    833                   iend = prodatqc%npvend(jobs,1)  
    834                   inum_obs = iend - ista + 1  
    835                   ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
    836                   DO iin=1,2   
    837                      DO ijn=1,2  
    838                                      
    839                                      
    840                         IF ( k1dint == 1 ) THEN  
    841                            CALL obs_int_z1d_spl( kpk, &  
    842                               &    zintt(iin,ijn,:,iobs),&  
    843                               &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    844                               &    zmask(iin,ijn,:,iobs))  
    845    
    846                         ENDIF  
    847         
    848                         CALL obs_level_search(kpk, &  
    849                             &        zgdept(iin,ijn,:,iobs),&  
    850                             &        inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
    851                             &         iv_indic)  
    852                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
    853                             &          prodatqc%var(1)%vdep(ista:iend),     &  
    854                             &          zintt(iin,ijn,:,iobs),            &  
    855                             &          zobs2k,interp_corner(iin,ijn,:), &  
    856                             &          zgdept(iin,ijn,:,iobs),         &  
    857                             &          zmask(iin,ijn,:,iobs) )       
    858           
    859                      ENDDO  
    860                   ENDDO  
    861               
    862                ENDIF  
    863         
    864                !-------------------------------------------------------------  
    865                ! Compute the horizontal interpolation for every profile level  
    866                !-------------------------------------------------------------  
    867               
    868                DO ikn=1,inum_obs  
    869                   iend=ista+ikn-1  
    870  
    871                   l_zweig(:,:,1) = 0._wp  
    872  
    873                   ! This code forces the horizontal weights to be   
    874                   ! zero IF the observation is below the bottom of the   
    875                   ! corners of the interpolation nodes, Or if it is in   
    876                   ! the mask. This is important for observations are near   
    877                   ! steep bathymetry  
    878                   DO iin=1,2  
    879                      DO ijn=1,2  
    880       
    881                         depth_loop1: DO ik=kpk,2,-1  
    882                            IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
    883                              
    884                               l_zweig(iin,ijn,1) = &   
    885                                  & zweig(iin,ijn,1) * &  
    886                                  & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
    887                                  &  - prodatqc%var(1)%vdep(iend)),0._wp)  
    888                              
    889                               EXIT depth_loop1  
    890                            ENDIF  
    891                         ENDDO depth_loop1  
    892       
    893                      ENDDO  
    894                   ENDDO  
    895     
    896                   CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &  
    897                   &          prodatqc%var(1)%vmod(iend:iend) )  
    898   
    899                ENDDO  
    900   
    901   
    902                DEALLOCATE(interp_corner,iv_indic)  
    903            
    904             ENDIF  
    905         
    906   
    907             ! Salinity   
    908            
    909             IF ( prodatqc%npvend(jobs,2) > 0 ) THEN   
    910      
    911                zobsk(:) = obfillflt  
    912      
    913                IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN  
    914      
    915                   IF ( idayend == 0 )  THEN  
    916                     
    917                      ! Daily averaged moored buoy (MRB) data  
    918                     
    919                      ! vertically interpolate all 4 corners  
    920                      ista = prodatqc%npvsta(iobs,2)  
    921                      iend = prodatqc%npvend(iobs,2)  
    922                      inum_obs = iend - ista + 1  
    923                      ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
    924        
    925                      DO iin=1,2  
    926                         DO ijn=1,2  
    927                                         
    928                                         
    929             
    930                            IF ( k1dint == 1 ) THEN  
    931                               CALL obs_int_z1d_spl( kpk, &  
    932                                  &     zinms(iin,ijn,:,iobs), &  
    933                                  &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    934                                  &     zmask(iin,ijn,:,iobs))  
    935                            ENDIF  
    936         
    937                            CALL obs_level_search(kpk, &  
    938                               &    zgdept(iin,ijn,:,iobs), &  
    939                               &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    940                               &    iv_indic)  
    941                            CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
    942                               &    prodatqc%var(2)%vdep(ista:iend), &  
    943                               &    zinms(iin,ijn,:,iobs), &  
    944                               &    zobs2k, interp_corner(iin,ijn,:), &  
    945                               &    zgdept(iin,ijn,:,iobs), &  
    946                               &    zmask(iin,ijn,:,iobs))  
    947         
    948                         ENDDO  
    949                      ENDDO  
    950                     
    951                     
    952                   ELSE  
    953                  
    954                      CALL ctl_stop( ' A nonzero' //     &  
    955                         &           ' number of profile T BUOY data should' // &  
    956                         &           ' only occur at the end of a given day' )  
    957      
    958                   ENDIF  
    959           
    960                ELSE   
    961                  
    962                   ! Point data  
    963       
    964                   ! vertically interpolate all 4 corners  
    965                   ista = prodatqc%npvsta(jobs,2)  
    966                   iend = prodatqc%npvend(jobs,2)  
    967                   inum_obs = iend - ista + 1  
    968                   ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
    969                     
    970                   DO iin=1,2      
    971                      DO ijn=1,2   
    972                                   
    973                                   
    974                         IF ( k1dint == 1 ) THEN  
    975                            CALL obs_int_z1d_spl( kpk, &  
    976                               &    zints(iin,ijn,:,iobs),&  
    977                               &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    978                               &    zmask(iin,ijn,:,iobs))  
    979    
    980                         ENDIF  
    981         
    982                         CALL obs_level_search(kpk, &  
    983                            &        zgdept(iin,ijn,:,iobs),&  
    984                            &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    985                            &         iv_indic)  
    986                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,  &  
    987                            &          prodatqc%var(2)%vdep(ista:iend),     &  
    988                            &          zints(iin,ijn,:,iobs),               &  
    989                            &          zobs2k,interp_corner(iin,ijn,:),     &  
    990                            &          zgdept(iin,ijn,:,iobs),              &  
    991                            &          zmask(iin,ijn,:,iobs) )       
    992           
    993                      ENDDO  
    994                   ENDDO  
    995               
    996                ENDIF  
    997         
    998                !-------------------------------------------------------------  
    999                ! Compute the horizontal interpolation for every profile level  
    1000                !-------------------------------------------------------------  
    1001               
    1002                DO ikn=1,inum_obs  
    1003                   iend=ista+ikn-1  
    1004  
    1005                   l_zweig(:,:,1) = 0._wp 
    1006     
    1007                   ! This code forces the horizontal weights to be   
    1008                   ! zero IF the observation is below the bottom of the   
    1009                   ! corners of the interpolation nodes, Or if it is in   
    1010                   ! the mask. This is important for observations are near   
    1011                   ! steep bathymetry  
    1012                   DO iin=1,2  
    1013                      DO ijn=1,2  
    1014       
    1015                         depth_loop2: DO ik=kpk,2,-1  
    1016                            IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
    1017                              
    1018                               l_zweig(iin,ijn,1) = &   
    1019                                  &  zweig(iin,ijn,1) * &  
    1020                                  &  MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
    1021                                  &  - prodatqc%var(2)%vdep(iend)),0._wp)  
    1022                              
    1023                               EXIT depth_loop2  
    1024                            ENDIF  
    1025                         ENDDO depth_loop2  
    1026       
    1027                      ENDDO  
    1028                   ENDDO  
    1029     
    1030                   CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &  
    1031                   &          prodatqc%var(2)%vmod(iend:iend) )  
    1032   
    1033                ENDDO  
    1034   
    1035   
    1036                DEALLOCATE(interp_corner,iv_indic)  
    1037            
    1038             ENDIF  
    1039            
    1040          ENDIF  
    1041         
    1042       END DO  
    1043       
    1044       ! Deallocate the data for interpolation  
    1045       DEALLOCATE( &  
    1046          & igrdi, &  
    1047          & igrdj, &  
    1048          & zglam, &  
    1049          & zgphi, &  
    1050          & zmask, &  
    1051          & zintt, &  
    1052          & zints, &  
    1053          & zgdept,& 
    1054          & zgdepw & 
    1055          & )  
    1056       ! At the end of the day also get interpolated means  
    1057       IF ( idayend == 0 ) THEN  
    1058          DEALLOCATE( &  
    1059             & zinmt,  &  
    1060             & zinms   &  
    1061             & )  
    1062       ENDIF  
    1063      
    1064       prodatqc%nprofup = prodatqc%nprofup + ipro   
    1065         
    1066    END SUBROUTINE obs_pro_sco_opt  
    1067   
    1068    SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,         & 
    1069       &                    kit000, kdaystp, psurf, psurfmask, & 
    1070       &                    k2dint, ldnightav ) 
     671   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            & 
     672      &                     kit000, kdaystp, psurf, psurfmask,   & 
     673      &                     k2dint, ldnightav, plamscl, pphiscl, & 
     674      &                     lindegrees ) 
    1071675 
    1072676      !!----------------------------------------------------------------------- 
     
    1090694      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    1091695      !! 
     696      !!    Two horizontal averaging schemes are also available: 
     697      !!        - weighted radial footprint        (k2dint = 5) 
     698      !!        - weighted rectangular footprint   (k2dint = 6) 
     699      !! 
    1092700      !! 
    1093701      !! ** Action  : 
     
    1096704      !!      ! 07-03 (A. Weaver) 
    1097705      !!      ! 15-02 (M. Martin) Combined routine for surface types 
     706      !!      ! 17-03 (M. Martin) Added horizontal averaging options 
    1098707      !!----------------------------------------------------------------------- 
    1099708 
     
    1117726         & psurfmask                   ! Land-sea mask 
    1118727      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 
     728      REAL(KIND=wp), INTENT(IN) :: & 
     729         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions 
     730         & pphiscl                     ! This is the full width (rather than half-width) 
     731      LOGICAL, INTENT(IN) :: & 
     732         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 
    1119733 
    1120734      !! * Local declarations 
     
    1125739      INTEGER :: isurf 
    1126740      INTEGER :: iobs 
     741      INTEGER :: imaxifp, imaxjfp 
     742      INTEGER :: imodi, imodj 
    1127743      INTEGER :: idayend 
    1128744      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    1129          & igrdi, & 
    1130          & igrdj 
     745         & igrdi,   & 
     746         & igrdj,   & 
     747         & igrdip1, & 
     748         & igrdjp1 
    1131749      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    1132750         & icount_night,      & 
     
    1136754      REAL(wp), DIMENSION(1) :: zext, zobsmask 
    1137755      REAL(wp) :: zdaystp 
    1138       REAL(wp), DIMENSION(2,2,1) :: & 
    1139          & zweig 
    1140756      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     757         & zweig,  & 
    1141758         & zmask,  & 
    1142759         & zsurf,  & 
    1143760         & zsurfm, & 
     761         & zsurftmp, & 
    1144762         & zglam,  & 
    1145          & zgphi 
     763         & zgphi,  & 
     764         & zglamf, & 
     765         & zgphif 
     766 
    1146767      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    1147768         & zintmp,  & 
     
    1155776      inrc = kt - kit000 + 2 
    1156777      isurf = surfdataqc%nsstp(inrc) 
     778 
     779      ! Work out the maximum footprint size for the  
     780      ! interpolation/averaging in model grid-points - has to be even. 
     781 
     782      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 
     783 
    1157784 
    1158785      IF ( ldnightav ) THEN 
     
    1221848 
    1222849      ALLOCATE( & 
    1223          & igrdi(2,2,isurf), & 
    1224          & igrdj(2,2,isurf), & 
    1225          & zglam(2,2,isurf), & 
    1226          & zgphi(2,2,isurf), & 
    1227          & zmask(2,2,isurf), & 
    1228          & zsurf(2,2,isurf)  & 
     850         & zweig(imaxifp,imaxjfp,1),      & 
     851         & igrdi(imaxifp,imaxjfp,isurf), & 
     852         & igrdj(imaxifp,imaxjfp,isurf), & 
     853         & zglam(imaxifp,imaxjfp,isurf), & 
     854         & zgphi(imaxifp,imaxjfp,isurf), & 
     855         & zmask(imaxifp,imaxjfp,isurf), & 
     856         & zsurf(imaxifp,imaxjfp,isurf), & 
     857         & zsurftmp(imaxifp,imaxjfp,isurf),  & 
     858         & zglamf(imaxifp+1,imaxjfp+1,isurf), & 
     859         & zgphif(imaxifp+1,imaxjfp+1,isurf), & 
     860         & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 
     861         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 
    1229862         & ) 
    1230863 
    1231864      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
    1232865         iobs = jobs - surfdataqc%nsurfup 
    1233          igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1 
    1234          igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1 
    1235          igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1 
    1236          igrdj(1,2,iobs) = surfdataqc%mj(jobs) 
    1237          igrdi(2,1,iobs) = surfdataqc%mi(jobs) 
    1238          igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1 
    1239          igrdi(2,2,iobs) = surfdataqc%mi(jobs) 
    1240          igrdj(2,2,iobs) = surfdataqc%mj(jobs) 
     866         DO ji = 0, imaxifp 
     867            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 
     868             
     869            !Deal with wrap around in longitude 
     870            IF ( imodi < 1      ) imodi = imodi + jpiglo 
     871            IF ( imodi > jpiglo ) imodi = imodi - jpiglo 
     872             
     873            DO jj = 0, imaxjfp 
     874               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 
     875               !If model values are out of the domain to the north/south then 
     876               !set them to be the edge of the domain 
     877               IF ( imodj < 1      ) imodj = 1 
     878               IF ( imodj > jpjglo ) imodj = jpjglo 
     879 
     880               igrdip1(ji+1,jj+1,iobs) = imodi 
     881               igrdjp1(ji+1,jj+1,iobs) = imodj 
     882                
     883               IF ( ji >= 1 .AND. jj >= 1 ) THEN 
     884                  igrdi(ji,jj,iobs) = imodi 
     885                  igrdj(ji,jj,iobs) = imodj 
     886               ENDIF 
     887                
     888            END DO 
     889         END DO 
    1241890      END DO 
    1242891 
    1243       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     892      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1244893         &                  igrdi, igrdj, glamt, zglam ) 
    1245       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     894      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1246895         &                  igrdi, igrdj, gphit, zgphi ) 
    1247       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     896      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1248897         &                  igrdi, igrdj, psurfmask, zmask ) 
    1249       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     898      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1250899         &                  igrdi, igrdj, psurf, zsurf ) 
     900      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     901         &                  igrdip1, igrdjp1, glamf, zglamf ) 
     902      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     903         &                  igrdip1, igrdjp1, gphif, zgphif ) 
    1251904 
    1252905      ! At the end of the day get interpolated means 
    1253       IF (ldnightav ) THEN 
    1254          IF ( idayend == 0 ) THEN 
    1255  
    1256             ALLOCATE( & 
    1257                & zsurfm(2,2,isurf)  & 
    1258                & ) 
    1259  
    1260             CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, & 
    1261                &               surfdataqc%vdmean(:,:), zsurfm ) 
    1262  
    1263          ENDIF 
     906      IF ( idayend == 0 .AND. ldnightav ) THEN 
     907 
     908         ALLOCATE( & 
     909            & zsurfm(imaxifp,imaxjfp,isurf)  & 
     910            & ) 
     911 
     912         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 
     913            &               surfdataqc%vdmean(:,:), zsurfm ) 
     914 
    1264915      ENDIF 
    1265916 
     
    1290941         zphi = surfdataqc%rphi(jobs) 
    1291942 
    1292          ! Get weights to interpolate the model value to the observation point 
    1293          CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    1294             &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    1295             &                   zmask(:,:,iobs), zweig, zobsmask ) 
    1296  
    1297          ! Interpolate the model field to the observation point 
    1298943         IF ( ldnightav .AND. idayend == 0 ) THEN 
    1299944            ! Night-time averaged data 
    1300             CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext ) 
     945            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 
    1301946         ELSE 
    1302             CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs),  zext ) 
     947            zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 
     948         ENDIF 
     949 
     950         IF ( k2dint <= 4 ) THEN 
     951 
     952            ! Get weights to interpolate the model value to the observation point 
     953            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     954               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     955               &                   zmask(:,:,iobs), zweig, zobsmask ) 
     956 
     957            ! Interpolate the model value to the observation point  
     958            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
     959 
     960         ELSE 
     961 
     962            ! Get weights to average the model SLA to the observation footprint 
     963            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
     964               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     965               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
     966               &                   zmask(:,:,iobs), plamscl, pphiscl, & 
     967               &                   lindegrees, zweig, zobsmask ) 
     968 
     969            ! Average the model SST to the observation footprint 
     970            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
     971               &              zweig, zsurftmp(:,:,iobs),  zext ) 
     972 
    1303973         ENDIF 
    1304974 
     
    1310980            surfdataqc%rmod(jobs,1) = zext(1) 
    1311981         ENDIF 
     982          
     983         IF ( zext(1) == obfillflt ) THEN 
     984            ! If the observation value is a fill value, set QC flag to bad 
     985            surfdataqc%nqc(jobs) = 4 
     986         ENDIF 
    1312987 
    1313988      END DO 
     
    1315990      ! Deallocate the data for interpolation 
    1316991      DEALLOCATE( & 
     992         & zweig, & 
    1317993         & igrdi, & 
    1318994         & igrdj, & 
     
    1320996         & zgphi, & 
    1321997         & zmask, & 
    1322          & zsurf  & 
     998         & zsurf, & 
     999         & zsurftmp, & 
     1000         & zglamf, & 
     1001         & zgphif, & 
     1002         & igrdip1,& 
     1003         & igrdjp1 & 
    13231004         & ) 
    13241005 
    13251006      ! At the end of the day also deallocate night-time mean array 
    1326       IF ( ldnightav ) THEN 
    1327          IF ( idayend == 0 ) THEN 
    1328             DEALLOCATE( & 
    1329                & zsurfm  & 
    1330                & ) 
    1331          ENDIF 
     1007      IF ( idayend == 0 .AND. ldnightav ) THEN 
     1008         DEALLOCATE( & 
     1009            & zsurfm  & 
     1010            & ) 
    13321011      ENDIF 
    13331012 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    r7646 r9023  
    2323   USE obs_oper           ! Observation operators 
    2424   USE lib_mpp, ONLY :   ctl_warn, ctl_stop 
     25   USE bdy_oce, ONLY : &        ! Boundary information 
     26      idx_bdy, nb_bdy, ln_bdy 
    2527 
    2628   IMPLICIT NONE 
     
    4042CONTAINS 
    4143 
    42    SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 
     44   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 
     45                            kqc_cutoff ) 
    4346      !!---------------------------------------------------------------------- 
    4447      !!                    ***  ROUTINE obs_pre_sla  *** 
     
    5760      !!        !  2015-02  (M. Martin) Combined routine for surface types. 
    5861      !!---------------------------------------------------------------------- 
     62      !! * Modules used 
    5963      USE par_oce             ! Ocean parameters 
    6064      USE dom_oce, ONLY       :   glamt, gphit, tmask, nproc   ! Geographical information 
     
    6367      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc   ! Subset of surface data not failing screening 
    6468      LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land 
    65       ! 
     69      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary 
     70      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
     71      !! * Local declarations 
     72      INTEGER :: iqc_cutoff = 255   ! cut off for QC value 
    6673      INTEGER :: iyea0        ! Initial date 
    6774      INTEGER :: imon0        !  - (year, month, day, hour, minute) 
     
    7683      INTEGER :: inlasobs     !  - close to land 
    7784      INTEGER :: igrdobs      !  - fail the grid search 
     85      INTEGER :: ibdysobs     !  - close to open boundary 
    7886                              ! Global counters for observations that 
    7987      INTEGER :: iotdobsmpp     !  - outside time domain 
     
    8290      INTEGER :: inlasobsmpp    !  - close to land 
    8391      INTEGER :: igrdobsmpp     !  - fail the grid search 
    84       LOGICAL, DIMENSION(:), ALLOCATABLE ::   llvalid            ! SLA data selection 
     92      INTEGER :: ibdysobsmpp  !  - close to open boundary 
     93      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
     94         & llvalid            ! SLA data selection 
    8595      INTEGER :: jobs         ! Obs. loop variable 
    8696      INTEGER :: jstp         ! Time loop variable 
     
    107117      ilansobs = 0 
    108118      inlasobs = 0 
     119      ibdysobs = 0  
     120 
     121      ! Set QC cutoff to optional value if provided 
     122      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 
    109123 
    110124      ! ----------------------------------------------------------------------- 
     
    140154         &                 tmask(:,:,1), surfdata%nqc,  & 
    141155         &                 iosdsobs,     ilansobs,     & 
    142          &                 inlasobs,     ld_nea        ) 
     156         &                 inlasobs,     ld_nea,       & 
     157         &                 ibdysobs,     ld_bound_reject, & 
     158         &                 iqc_cutoff                     ) 
    143159 
    144160      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
    145161      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
    146162      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
     163      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
    147164 
    148165      ! ----------------------------------------------------------------------- 
     
    155172      ALLOCATE( llvalid(surfdata%nsurf) ) 
    156173       
    157       ! We want all data which has qc flags <= 10 
    158  
    159       llvalid(:)  = ( surfdata%nqc(:)  <= 10 ) 
     174      ! We want all data which has qc flags <= iqc_cutoff 
     175 
     176      llvalid(:)  = ( surfdata%nqc(:)  <= iqc_cutoff ) 
    160177 
    161178      ! The actual copying 
     
    190207               &            inlasobsmpp 
    191208         ENDIF 
     209         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 
     210            &            ibdysobsmpp   
    192211         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', & 
    193212            &            surfdataqc%nsurfmpp 
     
    225244      &                     kpi, kpj, kpk, & 
    226245      &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  & 
    227       &                     ld_nea, kdailyavtypes ) 
     246      &                     ld_nea, ld_bound_reject, kdailyavtypes,  kqc_cutoff ) 
    228247 
    229248!!---------------------------------------------------------------------- 
     
    241260      !! 
    242261      !!---------------------------------------------------------------------- 
    243       USE par_oce           ! Ocean parameters 
    244       USE dom_oce, ONLY :   gdept_1d, nproc   ! Geographical information 
     262      !! * Modules used 
     263      USE par_oce             ! Ocean parameters 
     264      USE dom_oce, ONLY : &   ! Geographical information 
     265         & gdept_1d,             & 
     266         & nproc 
    245267 
    246268      !! * Arguments 
     
    250272      LOGICAL, INTENT(IN) :: ld_var2 
    251273      LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land 
     274      LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary 
    252275      INTEGER, INTENT(IN) :: kpi, kpj, kpk        ! Local domain sizes 
    253276      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
     
    261284         & pgphi1, & 
    262285         & pgphi2 
     286      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
    263287 
    264288      !! * Local declarations 
     289      INTEGER :: iqc_cutoff = 255   ! cut off for QC value 
    265290      INTEGER :: iyea0        ! Initial date 
    266291      INTEGER :: imon0        !  - (year, month, day, hour, minute) 
     
    277302      INTEGER :: inlav1obs    !  - close to land (variable 1) 
    278303      INTEGER :: inlav2obs    !  - close to land (variable 2) 
     304      INTEGER :: ibdyv1obs    !  - boundary (variable 1)  
     305      INTEGER :: ibdyv2obs    !  - boundary (variable 2)       
    279306      INTEGER :: igrdobs      !  - fail the grid search 
    280307      INTEGER :: iuvchku      !  - reject u if v rejected and vice versa 
     
    288315      INTEGER :: inlav1obsmpp !  - close to land (variable 1) 
    289316      INTEGER :: inlav2obsmpp !  - close to land (variable 2) 
     317      INTEGER :: ibdyv1obsmpp !  - boundary (variable 1)  
     318      INTEGER :: ibdyv2obsmpp !  - boundary (variable 2)       
    290319      INTEGER :: igrdobsmpp   !  - fail the grid search 
    291320      INTEGER :: iuvchkumpp   !  - reject var1 if var2 rejected and vice versa 
     
    322351      inlav1obs = 0 
    323352      inlav2obs = 0 
    324       iuvchku  = 0 
    325       iuvchkv = 0 
     353      ibdyv1obs = 0 
     354      ibdyv2obs = 0 
     355      iuvchku   = 0 
     356      iuvchkv   = 0 
     357 
     358 
     359      ! Set QC cutoff to optional value if provided 
     360      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 
    326361 
    327362      ! ----------------------------------------------------------------------- 
     
    335370            &              profdata%nday,    profdata%nhou, profdata%nmin, & 
    336371            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    337             &              iotdobs, kdailyavtypes = kdailyavtypes ) 
     372            &              iotdobs, kdailyavtypes = kdailyavtypes,         & 
     373            &              kqc_cutoff = iqc_cutoff ) 
    338374      ELSE 
    339375         CALL obs_coo_tim_prof( icycle, & 
     
    342378            &              profdata%nday,    profdata%nhou, profdata%nmin, & 
    343379            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    344             &              iotdobs ) 
     380            &              iotdobs,          kqc_cutoff = iqc_cutoff ) 
    345381      ENDIF 
    346382 
     
    359395 
    360396      ! ----------------------------------------------------------------------- 
    361       ! Reject all observations for profiles with nqc > 10 
    362       ! ----------------------------------------------------------------------- 
    363  
    364       CALL obs_pro_rej( profdata ) 
     397      ! Reject all observations for profiles with nqc > iqc_cutoff 
     398      ! ----------------------------------------------------------------------- 
     399 
     400      CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff ) 
    365401 
    366402      ! ----------------------------------------------------------------------- 
     
    381417         &                 gdept_1d,              zmask1,               & 
    382418         &                 profdata%nqc,          profdata%var(1)%nvqc, & 
    383          &                 iosdv1obs,              ilanv1obs,           & 
    384          &                 inlav1obs,              ld_nea                ) 
     419         &                 iosdv1obs,             ilanv1obs,            & 
     420         &                 inlav1obs,             ld_nea,               & 
     421         &                 ibdyv1obs,             ld_bound_reject,      & 
     422         &                 iqc_cutoff       ) 
    385423 
    386424      CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 
    387425      CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 
    388426      CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 
     427      CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp ) 
    389428 
    390429      ! Variable 2 
     
    400439         &                 gdept_1d,              zmask2,               & 
    401440         &                 profdata%nqc,          profdata%var(2)%nvqc, & 
    402          &                 iosdv2obs,              ilanv2obs,           & 
    403          &                 inlav2obs,              ld_nea                ) 
     441         &                 iosdv2obs,             ilanv2obs,            & 
     442         &                 inlav2obs,             ld_nea,               & 
     443         &                 ibdyv2obs,             ld_bound_reject,      & 
     444         &                 iqc_cutoff       ) 
    404445 
    405446      CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 
    406447      CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 
    407448      CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 
     449      CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 
    408450 
    409451      ! ----------------------------------------------------------------------- 
     
    412454 
    413455      IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    414          CALL obs_uv_rej( profdata, iuvchku, iuvchkv ) 
     456         CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 
    415457         CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 
    416458         CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) 
     
    429471      END DO 
    430472 
    431       ! We want all data which has qc flags = 0 
    432  
    433       llvalid%luse(:) = ( profdata%nqc(:)  <= 10 ) 
     473      ! We want all data which has qc flags <= iqc_cutoff 
     474 
     475      llvalid%luse(:) = ( profdata%nqc(:)  <= iqc_cutoff ) 
    434476      DO jvar = 1,profdata%nvar 
    435          llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 ) 
     477         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff ) 
    436478      END DO 
    437479 
     
    475517               &            iuvchku 
    476518         ENDIF 
     519         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',& 
     520               &            ibdyv1obsmpp 
    477521         WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', & 
    478522            &            prodatqc%nvprotmpp(1) 
     
    492536               &            iuvchkv 
    493537         ENDIF 
     538         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 
     539               &            ibdyv2obsmpp 
    494540         WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', & 
    495541            &            prodatqc%nvprotmpp(2) 
     
    644690            &        .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN 
    645691            kobsstp(jobs) = -1 
    646             kobsqc(jobs)  = kobsqc(jobs) + 11 
     692            kobsqc(jobs)  = IBSET(kobsqc(jobs),13) 
    647693            kotdobs       = kotdobs + 1 
    648694            CYCLE 
     
    695741         IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) & 
    696742            & .OR.( kobsstp(jobs) > nitend ) ) THEN 
    697             kobsqc(jobs) = kobsqc(jobs) + 12 
     743            kobsqc(jobs) = IBSET(kobsqc(jobs),13) 
    698744            kotdobs = kotdobs + 1 
    699745            CYCLE 
     
    739785      &                    kobsno,                                        & 
    740786      &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   & 
    741       &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes ) 
     787      &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes, & 
     788      &                    kqc_cutoff ) 
    742789      !!---------------------------------------------------------------------- 
    743790      !!                    ***  ROUTINE obs_coo_tim *** 
     
    783830      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    784831         & kdailyavtypes    ! Types for daily averages 
     832      INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff           ! QC cutoff value 
     833 
    785834      !! * Local declarations 
    786835      INTEGER :: jobs 
     836      INTEGER :: iqc_cutoff=255 
    787837 
    788838      !----------------------------------------------------------------------- 
     
    803853         DO jobs = 1, kobsno 
    804854             
    805             IF ( kobsqc(jobs) <= 10 ) THEN 
     855            IF ( kobsqc(jobs) <= iqc_cutoff ) THEN 
    806856                
    807857               IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.& 
    808858                  & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN 
    809                   kobsqc(jobs) = kobsqc(jobs) + 14 
     859                  kobsqc(jobs) = IBSET(kobsqc(jobs),13) 
    810860                  kotdobs      = kotdobs + 1 
    811861                  CYCLE 
     
    850900      DO jobs = 1, kobsno 
    851901         IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN 
    852             kobsqc(jobs) = kobsqc(jobs) + 18 
     902            kobsqc(jobs) = IBSET(kobsqc(jobs),12) 
    853903            kgrdobs = kgrdobs + 1 
    854904         ENDIF 
     
    861911      &                       plam,   pphi,    pmask,            & 
    862912      &                       kobsqc, kosdobs, klanobs,          & 
    863       &                       knlaobs,ld_nea                     ) 
     913      &                       knlaobs,ld_nea,                    & 
     914      &                       kbdyobs,ld_bound_reject,           & 
     915      &                       kqc_cutoff                         ) 
    864916      !!---------------------------------------------------------------------- 
    865917      !!                    ***  ROUTINE obs_coo_spc_2d  *** 
     
    894946      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 
    895947         & kobsqc             ! Observation quality control 
    896       INTEGER, INTENT(INOUT) :: kosdobs   ! Observations outside space domain 
    897       INTEGER, INTENT(INOUT) :: klanobs   ! Observations within a model land cell 
    898       INTEGER, INTENT(INOUT) :: knlaobs   ! Observations near land 
    899       LOGICAL, INTENT(IN) :: ld_nea       ! Flag observations near land 
     948      INTEGER, INTENT(INOUT) :: kosdobs          ! Observations outside space domain 
     949      INTEGER, INTENT(INOUT) :: klanobs          ! Observations within a model land cell 
     950      INTEGER, INTENT(INOUT) :: knlaobs          ! Observations near land 
     951      INTEGER, INTENT(INOUT) :: kbdyobs          ! Observations near boundary 
     952      LOGICAL, INTENT(IN)    :: ld_nea           ! Flag observations near land 
     953      LOGICAL, INTENT(IN)    :: ld_bound_reject  ! Flag observations near open boundary  
     954      INTEGER, INTENT(IN)    :: kqc_cutoff       ! Cutoff QC value 
     955 
    900956      !! * Local declarations 
    901957      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    902958         & zgmsk              ! Grid mask 
     959 
     960      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
     961         & zbmsk              ! Boundary mask 
     962      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
    903963      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    904964         & zglam, &           ! Model longitude at grid points 
     
    917977         ! For invalid points use 2,2 
    918978 
    919          IF ( kobsqc(jobs) >= 10 ) THEN 
     979         IF ( kobsqc(jobs) >= kqc_cutoff ) THEN 
    920980 
    921981            igrdi(1,1,jobs) = 1 
     
    9421002 
    9431003      END DO 
     1004 
     1005      IF (ln_bdy) THEN 
     1006        ! Create a mask grid points in boundary rim 
     1007        IF (ld_bound_reject) THEN 
     1008           zbdymask(:,:) = 1.0_wp 
     1009           DO ji = 1, nb_bdy 
     1010              DO jj = 1, idx_bdy(ji)%nblen(1) 
     1011                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     1012              ENDDO 
     1013           ENDDO 
     1014 
     1015           CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 
     1016        ENDIF 
     1017      ENDIF 
     1018 
    9441019       
    9451020      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) 
     
    9501025 
    9511026         ! Skip bad observations 
    952          IF ( kobsqc(jobs) >= 10 ) CYCLE 
     1027         IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE 
    9531028 
    9541029         ! Flag if the observation falls outside the model spatial domain 
     
    9571032            &  .OR. ( pobsphi(jobs) <  -90. ) & 
    9581033            &  .OR. ( pobsphi(jobs) >   90. ) ) THEN 
    959             kobsqc(jobs) = kobsqc(jobs) + 11 
     1034            kobsqc(jobs) = IBSET(kobsqc(jobs),11) 
    9601035            kosdobs = kosdobs + 1 
    9611036            CYCLE 
     
    9641039         ! Flag if the observation falls with a model land cell 
    9651040         IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
    966             kobsqc(jobs) = kobsqc(jobs)  + 12 
     1041            kobsqc(jobs) = IBSET(kobsqc(jobs),10) 
    9671042            klanobs = klanobs + 1 
    9681043            CYCLE 
     
    9781053               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 
    9791054                  & .AND. & 
    980                   & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) & 
    981                   & ) THEN 
     1055                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) & 
     1056                  & < 1.0e-6_wp ) ) THEN 
    9821057                  lgridobs = .TRUE. 
    9831058                  iig = ji 
     
    9921067         IF (lgridobs) THEN 
    9931068            IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
    994                kobsqc(jobs) = kobsqc(jobs) + 12 
     1069               kobsqc(jobs) = IBSET(kobsqc(jobs),10) 
    9951070               klanobs = klanobs + 1 
    9961071               CYCLE 
     
    10001075         ! Flag if the observation falls is close to land 
    10011076         IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN 
    1002             IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 14 
    10031077            knlaobs = knlaobs + 1 
    1004             CYCLE 
     1078            IF (ld_nea) THEN 
     1079               kobsqc(jobs) = IBSET(kobsqc(jobs),9) 
     1080               CYCLE 
     1081            ENDIF 
     1082         ENDIF 
     1083 
     1084         IF (ln_bdy) THEN 
     1085         ! Flag if the observation falls close to the boundary rim 
     1086           IF (ld_bound_reject) THEN 
     1087              IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
     1088                 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 
     1089                 kbdyobs = kbdyobs + 1 
     1090                 CYCLE 
     1091              ENDIF 
     1092              ! for observations on the grid... 
     1093              IF (lgridobs) THEN 
     1094                 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
     1095                    kobsqc(jobs) = IBSET(kobsqc(jobs),8) 
     1096                    kbdyobs = kbdyobs + 1 
     1097                    CYCLE 
     1098                 ENDIF 
     1099              ENDIF 
     1100            ENDIF 
    10051101         ENDIF 
    10061102             
     
    10151111      &                       plam,    pphi,    pdep,    pmask, & 
    10161112      &                       kpobsqc, kobsqc,  kosdobs,        & 
    1017       &                       klanobs, knlaobs, ld_nea          ) 
     1113      &                       klanobs, knlaobs, ld_nea,         & 
     1114      &                       kbdyobs, ld_bound_reject,         & 
     1115      &                       kqc_cutoff                        ) 
    10181116      !!---------------------------------------------------------------------- 
    10191117      !!                    ***  ROUTINE obs_coo_spc_3d  *** 
     
    10771175      INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell 
    10781176      INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land 
     1177      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary 
    10791178      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land 
     1179      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary 
     1180      INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value 
     1181 
    10801182      !! * Local declarations 
    10811183      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
    10821184         & zgmsk              ! Grid mask 
     1185      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 
     1186         & zbmsk              ! Boundary mask 
     1187      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
    10831188      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
    10841189         & zgdepw 
     
    11001205         ! For invalid points use 2,2 
    11011206 
    1102          IF ( kpobsqc(jobs) >= 10 ) THEN 
     1207         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN 
    11031208             
    11041209            igrdi(1,1,jobs) = 1 
     
    11251230          
    11261231      END DO 
     1232 
     1233      IF (ln_bdy) THEN 
     1234        ! Create a mask grid points in boundary rim 
     1235        IF (ld_bound_reject) THEN            
     1236           zbdymask(:,:) = 1.0_wp 
     1237           DO ji = 1, nb_bdy 
     1238              DO jj = 1, idx_bdy(ji)%nblen(1) 
     1239                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     1240              ENDDO 
     1241           ENDDO 
     1242        ENDIF 
     1243    
     1244        CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 
     1245      ENDIF 
    11271246       
    11281247      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 
    11291248      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 
    11301249      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 
    1131       IF ( .NOT.( ln_zps .OR. ln_zco ) ) THEN 
    1132         ! Need to know the bathy depth for each observation for sco 
    1133         CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 
     1250      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 
    11341251        &                     zgdepw ) 
    1135       ENDIF