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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5120 for trunk/NEMOGCM/NEMO/OPA_SRC – NEMO

Ignore:
Timestamp:
2015-03-03T17:11:55+01:00 (9 years ago)
Author:
acc
Message:

#1473. Re-organisation and optimisation of ice shelf cavity option. This commit merges changes from the dev_r5094_UKMO_ISFCLEAN branch onto the trunk. Results will change, even with ln_isfcav=F, due to a return to original definitions of the vertical metrics. All changes have been reviewed and SETTE tested.

Location:
trunk/NEMOGCM/NEMO/OPA_SRC
Files:
39 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r4998 r5120  
    746746 
    747747 
    748             IF( ln_zps .AND. .NOT. lk_c1d ) & 
    749                &  CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    750                &                rhd, gru , grv, aru, arv, gzu, gzv, ge3ru, ge3rv,  &             ! 
    751                &                gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     748            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      & 
     749               &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient 
     750               &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level 
     751            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      & 
     752               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,    &    ! Partial steps for top cell (ISF) 
     753               &                                  rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     754               &                           gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    752755 
    753756#if defined key_zdfkpp 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r5107 r5120  
    103103      END DO 
    104104      IF( .NOT.lk_vvl ) THEN 
    105          DO ji=1,jpi 
    106             DO jj=1,jpj 
    107                zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
    108             END DO 
    109          END DO 
     105         IF ( ln_isfcav ) THEN 
     106            DO ji=1,jpi 
     107               DO jj=1,jpj 
     108                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     109               END DO 
     110            END DO 
     111         ELSE 
     112            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     113         END IF 
    110114      END IF 
    111115      !                                          
    … …  
    125129      END DO 
    126130      IF( .NOT.lk_vvl ) THEN 
    127          DO ji=1,jpi 
    128             DO jj=1,jpj 
    129                zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
    130             END DO 
    131          END DO 
     131         IF ( ln_isfcav ) THEN 
     132            DO ji=1,jpi 
     133               DO jj=1,jpj 
     134                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     135               END DO 
     136            END DO 
     137         ELSE 
     138            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     139         END IF 
    132140      END IF 
    133141      !     
    … …  
    155163      END DO 
    156164      IF( .NOT.lk_vvl ) THEN 
    157          DO ji=1,jpi 
    158             DO jj=1,jpj 
    159                ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)  
    160                zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)  
    161             END DO 
    162          END DO 
     165         IF ( ln_isfcav ) THEN 
     166            DO ji=1,jpi 
     167               DO jj=1,jpj 
     168                  ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)  
     169                  zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)  
     170               END DO 
     171            END DO 
     172         ELSE 
     173            ztemp = ztemp + zarea_ssh(:,:) * tsn(:,:,1,jp_tem)  
     174            zsal  = zsal  + zarea_ssh(:,:) * tsn(:,:,1,jp_sal)  
     175         END IF 
    163176      ENDIF 
    164177      IF( lk_mpp ) THEN   

  • trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r4990 r5120  
    9696      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )                               ! heat fluxes 
    9797      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )                               ! salt fluxes 
    98       ! Add runoff heat & salt input 
     98      ! Add runoff    heat & salt input 
    9999      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * surf(:,:) ) 
    100100      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) ) 
    101       ! Add geothermal ice shelf 
     101      ! Add ice shelf heat & salt input 
    102102      IF( nn_isf .GE. 1 )  THEN 
    103103          z_frc_trd_t = z_frc_trd_t & 
    … …  
    112112      ! 
    113113      IF( .NOT. lk_vvl ) THEN 
    114          z2d0=0.0_wp ; z2d1=0.0_wp 
    115          DO ji=1,jpi 
    116             DO jj=1,jpj 
    117               z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem) 
    118               z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal) 
     114         IF ( ln_isfcav ) THEN 
     115            DO ji=1,jpi 
     116               DO jj=1,jpj 
     117                  z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem) 
     118                  z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal) 
     119               ENDDO 
    119120            ENDDO 
    120          ENDDO 
     121         ELSE 
     122            z2d0(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) 
     123            z2d1(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) 
     124         END IF 
    121125         z_wn_trd_t = - glob_sum( z2d0 )  
    122126         z_wn_trd_s = - glob_sum( z2d1 ) 
    … …  
    144148      ! heat & salt content variation (associated with ssh) 
    145149      IF( .NOT. lk_vvl ) THEN 
    146          z2d0 = 0._wp   ;   z2d1 = 0._wp 
    147          DO ji = 1, jpi 
    148             DO jj = 1, jpj 
    149               z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) )  
    150               z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) )  
     150         IF ( ln_isfcav ) THEN 
     151            DO ji = 1, jpi 
     152               DO jj = 1, jpj 
     153                  z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) )  
     154                  z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) )  
     155               END DO 
    151156            END DO 
    152          END DO 
     157         ELSE 
     158            z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) )  
     159            z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) )  
     160         END IF 
    153161         z_ssh_hc = glob_sum( z2d0 )  
    154162         z_ssh_sc = glob_sum( z2d1 )  
    … …  
    277285          frc_s = 0._wp                                           ! salt content   -    -   -    -         
    278286          IF( .NOT. lk_vvl ) THEN 
    279              DO ji=1,jpi 
    280                 DO jj=1,jpj 
    281                    ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh 
    282                    ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh 
     287             IF ( ln_isfcav ) THEN 
     288                DO ji=1,jpi 
     289                   DO jj=1,jpj 
     290                      ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh 
     291                      ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh 
     292                   ENDDO 
    283293                ENDDO 
    284              ENDDO 
     294             ELSE 
     295                ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh 
     296                ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh 
     297             END IF 
    285298             frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface 
    286299             frc_wn_s = 0._wp                                       ! initial salt content misfit due to free surface 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r4990 r5120  
    262262 
    263263   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
     264   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts 
    264265 
    265266   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4) 
    … …  
    332333   INTEGER FUNCTION dom_oce_alloc() 
    333334      !!---------------------------------------------------------------------- 
    334       INTEGER, DIMENSION(11) :: ierr 
     335      INTEGER, DIMENSION(12) :: ierr 
    335336      !!---------------------------------------------------------------------- 
    336337      ierr(:) = 0 
    … …  
    400401         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(11) ) 
    401402 
     403      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 
     404 
    402405#if defined key_noslip_accurate 
    403       ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(11) ) 
     406      ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) ) 
    404407#endif 
    405408      ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r4990 r5120  
    281281      CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
    282282 
     283      ! 3. Ocean/land mask at wu-, wv- and w points  
     284      !---------------------------------------------- 
     285      wmask (:,:,1) = tmask(:,:,1) ! ???????? 
     286      wumask(:,:,1) = umask(:,:,1) ! ???????? 
     287      wvmask(:,:,1) = vmask(:,:,1) ! ???????? 
     288      DO jk=2,jpk 
     289         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1) 
     290         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)    
     291         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1) 
     292      END DO 
    283293 
    284294      ! 4. ocean/land mask for the elliptic equation 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5107 r5120  
    88   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: 
    99   !!                                          vvl option includes z_star and z_tilde coordinates 
     10   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    1011   !!---------------------------------------------------------------------- 
    1112   !!   'key_vvl'                              variable volume 
    … …  
    125126      INTEGER ::   ji,jj,jk 
    126127      INTEGER ::   ii0, ii1, ij0, ij1 
     128      REAL(wp)::   zcoef 
    127129      !!---------------------------------------------------------------------- 
    128130      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init') 
    … …  
    164166      ! t- and w- points depth 
    165167      ! ---------------------- 
     168      ! set the isf depth as it is in the initial step 
    166169      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
    167170      fsdepw_n(:,:,1) = 0.0_wp 
    … …  
    169172      fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 
    170173      fsdepw_b(:,:,1) = 0.0_wp 
    171       DO jj = 1,jpj 
    172          DO ji = 1,jpi 
    173             DO jk = 2,mikt(ji,jj)-1 
    174                fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 
    175                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    176                fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 
    177                fsdept_b(ji,jj,jk) = gdept_0(ji,jj,jk) 
    178                fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    179             END DO 
    180             IF (mikt(ji,jj) .GT. 1) THEN 
    181                jk = mikt(ji,jj) 
    182                fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 
    183                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    184                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    185                fsdept_b(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_b(ji,jj,jk) 
    186                fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    187             END IF 
    188             DO jk = mikt(ji,jj)+1, jpk 
    189                fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 
     174 
     175      DO jk = 2, jpk 
     176         DO jj = 1,jpj 
     177            DO ji = 1,jpi 
     178              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     179                                                     ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     180                                                     ! 0.5 where jk = mikt   
     181               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    190182               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
    191                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    192                fsdept_b(ji,jj,jk) = fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk) 
     183               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  & 
     184                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk))  
     185               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 
    193186               fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 
     187               fsdept_b(ji,jj,jk) =      zcoef  * ( fsdepw_b(ji,jj,jk  ) + 0.5 * fse3w_b(ji,jj,jk))  & 
     188                   &                + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) +       fse3w_b(ji,jj,jk))  
    194189            END DO 
    195190         END DO 
    … …  
    589584      !! * Local declarations 
    590585      INTEGER                             :: ji,jj,jk       ! dummy loop indices 
     586      REAL(wp)                            :: zcoef 
    591587      !!---------------------------------------------------------------------- 
    592588 
    … …  
    635631      ! t- and w- points depth 
    636632      ! ---------------------- 
     633      ! set the isf depth as it is in the initial step 
    637634      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
    638635      fsdepw_n(:,:,1) = 0.0_wp 
    639636      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 
    640       DO jj = 1,jpj 
    641          DO ji = 1,jpi 
    642             DO jk = 2,mikt(ji,jj)-1 
    643                fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 
    644                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    645                fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 
    646             END DO 
    647             IF (mikt(ji,jj) .GT. 1) THEN 
    648                jk = mikt(ji,jj) 
    649                fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 
    650                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    651                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    652             END IF 
    653             DO jk = mikt(ji,jj)+1, jpk 
    654                fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 
     637 
     638      DO jk = 2, jpk 
     639         DO jj = 1,jpj 
     640            DO ji = 1,jpi 
     641              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     642                                                                 ! 1 for jk = mikt 
     643               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    655644               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
    656                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
     645               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  & 
     646                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk))  
     647               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 
    657648            END DO 
    658649         END DO 
    659650      END DO 
     651 
    660652      ! Local depth and Inverse of the local depth of the water column at u- and v- points 
    661653      ! ---------------------------------------------------------------------------------- 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5118 r5120  
    1717   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function 
    1818   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case   
     19   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye   
    1920   !!---------------------------------------------------------------------- 
    2021 
    … …  
    3536   USE oce               ! ocean variables 
    3637   USE dom_oce           ! ocean domain 
    37    USE sbc_oce           ! surface variable (isf) 
    3838   USE closea            ! closed seas 
    3939   USE c1d               ! 1D vertical configuration 
    … …  
    298298      ENDIF 
    299299 
     300      IF ( ln_isfcav ) THEN 
    300301! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 
    301302! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
    302       DO jk = 1, jpkm1 
    303          e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
    304       END DO 
    305       e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO 
    306  
    307       DO jk = 2, jpk 
    308          e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)  
    309       END DO 
    310       e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))  
     303         DO jk = 1, jpkm1 
     304            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
     305         END DO 
     306         e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO 
     307 
     308         DO jk = 2, jpk 
     309            e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)  
     310         END DO 
     311         e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))  
     312      END IF 
    311313 
    312314!!gm BUG in s-coordinate this does not work! 
    … …  
    472474         ! 
    473475         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 
    474          IF( cp_cfg == "isomip" ) THEN  
    475            !  
    476            risfdep(:,:)=200.e0  
    477            misfdep(:,:)=1  
    478            ij0 = 1 ; ij1 = 40  
    479            DO jj = mj0(ij0), mj1(ij1)  
    480               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp  
    481                 END DO  
     476         IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN  
     477            risfdep(:,:)=200.e0  
     478            misfdep(:,:)=1  
     479            ij0 = 1 ; ij1 = 40  
     480            DO jj = mj0(ij0), mj1(ij1)  
     481               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp  
     482            END DO  
    482483            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp  
    483            !  
    484          ELSEIF ( cp_cfg == "isomip2" ) THEN 
     484         !  
     485         ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN 
    485486         !  
    486487            risfdep(:,:)=0.e0 
    … …  
    540541            END IF 
    541542            CALL iom_close( inum ) 
    542             !   
     543            !                                                 
    543544            risfdep(:,:)=0._wp          
    544545            misfdep(:,:)=1              
    … …  
    588589      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==! 
    589590         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 
    590          WHERE (bathy == risfdep) 
    591             bathy   = 0.0_wp ; risfdep = 0.0_wp 
    592          END WHERE 
     591         IF ( ln_isfcav ) THEN 
     592            WHERE (bathy == risfdep) 
     593               bathy   = 0.0_wp ; risfdep = 0.0_wp 
     594            END WHERE 
     595         END IF 
    593596         ! end patch 
    594597         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level 
    … …  
    965968      !!---------------------------------------------------------------------- 
    966969      !! 
     970      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     971      INTEGER  ::   ik, it           ! temporary integers 
     972      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
     973      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     974      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
     975      REAL(wp) ::   zmax             ! Maximum depth 
     976      REAL(wp) ::   zdiff            ! temporary scalar 
     977      REAL(wp) ::   zrefdep          ! temporary scalar 
     978      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
     979      !!--------------------------------------------------------------------- 
     980      ! 
     981      IF( nn_timing == 1 )  CALL timing_start('zgr_zps') 
     982      ! 
     983      CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     984      ! 
     985      IF(lwp) WRITE(numout,*) 
     986      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps' 
     987      IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
     988      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
     989 
     990      ll_print = .FALSE.                   ! Local variable for debugging 
     991       
     992      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
     993         WRITE(numout,*) 
     994         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
     995         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
     996      ENDIF 
     997 
     998 
     999      ! bathymetry in level (from bathy_meter) 
     1000      ! =================== 
     1001      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 
     1002      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
     1003      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
     1004      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level 
     1005      END WHERE 
     1006 
     1007      ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
     1008      ! find the number of ocean levels such that the last level thickness 
     1009      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 
     1010      ! e3t_1d is the reference level thickness 
     1011      DO jk = jpkm1, 1, -1 
     1012         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1013         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
     1014      END DO 
     1015 
     1016      IF ( ln_isfcav ) CALL zgr_isf 
     1017 
     1018      ! Scale factors and depth at T- and W-points 
     1019      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
     1020         gdept_0(:,:,jk) = gdept_1d(jk) 
     1021         gdepw_0(:,:,jk) = gdepw_1d(jk) 
     1022         e3t_0  (:,:,jk) = e3t_1d  (jk) 
     1023         e3w_0  (:,:,jk) = e3w_1d  (jk) 
     1024      END DO 
     1025      !  
     1026      DO jj = 1, jpj 
     1027         DO ji = 1, jpi 
     1028            ik = mbathy(ji,jj) 
     1029            IF( ik > 0 ) THEN               ! ocean point only 
     1030               ! max ocean level case 
     1031               IF( ik == jpkm1 ) THEN 
     1032                  zdepwp = bathy(ji,jj) 
     1033                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
     1034                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
     1035                  e3t_0(ji,jj,ik  ) = ze3tp 
     1036                  e3t_0(ji,jj,ik+1) = ze3tp 
     1037                  e3w_0(ji,jj,ik  ) = ze3wp 
     1038                  e3w_0(ji,jj,ik+1) = ze3tp 
     1039                  gdepw_0(ji,jj,ik+1) = zdepwp 
     1040                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
     1041                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
     1042                  ! 
     1043               ELSE                         ! standard case 
     1044                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
     1045                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     1046                  ENDIF 
     1047!gm Bug?  check the gdepw_1d 
     1048                  !       ... on ik 
     1049                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     1050                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     1051                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     1052                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     1053                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     1054                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
     1055                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     1056                  !       ... on ik+1 
     1057                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1058                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1059                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
     1060               ENDIF 
     1061            ENDIF 
     1062         END DO 
     1063      END DO 
     1064      ! 
     1065      it = 0 
     1066      DO jj = 1, jpj 
     1067         DO ji = 1, jpi 
     1068            ik = mbathy(ji,jj) 
     1069            IF( ik > 0 ) THEN               ! ocean point only 
     1070               e3tp (ji,jj) = e3t_0(ji,jj,ik) 
     1071               e3wp (ji,jj) = e3w_0(ji,jj,ik) 
     1072               ! test 
     1073               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
     1074               IF( zdiff <= 0._wp .AND. lwp ) THEN  
     1075                  it = it + 1 
     1076                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
     1077                  WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
     1078                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
     1079                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
     1080               ENDIF 
     1081            ENDIF 
     1082         END DO 
     1083      END DO 
     1084      ! 
     1085      IF ( ln_isfcav ) THEN 
     1086      ! (ISF) Definition of e3t, u, v, w for ISF case 
     1087         DO jj = 1, jpj  
     1088            DO ji = 1, jpi  
     1089               ik = misfdep(ji,jj)  
     1090               IF( ik > 1 ) THEN               ! ice shelf point only  
     1091                  IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
     1092                  gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
     1093!gm Bug?  check the gdepw_0  
     1094               !       ... on ik  
     1095                  gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
     1096                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
     1097                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
     1098                  e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1099                  e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1100 
     1101                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
     1102                     e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1103                  ENDIF  
     1104               !       ... on ik / ik-1  
     1105                  e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1106                  e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
     1107! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
     1108                  gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1109               ENDIF  
     1110            END DO  
     1111         END DO  
     1112      !  
     1113         it = 0  
     1114         DO jj = 1, jpj  
     1115            DO ji = 1, jpi  
     1116               ik = misfdep(ji,jj)  
     1117               IF( ik > 1 ) THEN               ! ice shelf point only  
     1118                  e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
     1119                  e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
     1120               ! test  
     1121                  zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
     1122                  IF( zdiff <= 0. .AND. lwp ) THEN   
     1123                     it = it + 1  
     1124                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
     1125                     WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
     1126                     WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
     1127                     WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
     1128                  ENDIF  
     1129               ENDIF  
     1130            END DO  
     1131         END DO  
     1132      END IF 
     1133      ! END (ISF) 
     1134 
     1135      ! Scale factors and depth at U-, V-, UW and VW-points 
     1136      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1137         e3u_0 (:,:,jk) = e3t_1d(jk) 
     1138         e3v_0 (:,:,jk) = e3t_1d(jk) 
     1139         e3uw_0(:,:,jk) = e3w_1d(jk) 
     1140         e3vw_0(:,:,jk) = e3w_1d(jk) 
     1141      END DO 
     1142      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
     1143         DO jj = 1, jpjm1 
     1144            DO ji = 1, fs_jpim1   ! vector opt. 
     1145               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
     1146               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
     1147               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
     1148               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
     1149            END DO 
     1150         END DO 
     1151      END DO 
     1152      IF ( ln_isfcav ) THEN 
     1153      ! (ISF) define e3uw (adapted for 2 cells in the water column) 
     1154      ! Need to test if the modification of only mikt and mbkt levels is enough 
     1155         DO jk = 2,jpk                           
     1156            DO jj = 1, jpjm1  
     1157               DO ji = 1, fs_jpim1   ! vector opt.  
     1158                  e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) & 
     1159                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) ) 
     1160                  e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) & 
     1161                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) ) 
     1162               END DO  
     1163            END DO  
     1164         END DO 
     1165      END IF 
     1166       
     1167      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
     1168      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
     1169      ! 
     1170      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1171         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
     1172         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
     1173         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
     1174         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
     1175      END DO 
     1176       
     1177      ! Scale factor at F-point 
     1178      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1179         e3f_0(:,:,jk) = e3t_1d(jk) 
     1180      END DO 
     1181      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
     1182         DO jj = 1, jpjm1 
     1183            DO ji = 1, fs_jpim1   ! vector opt. 
     1184               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
     1185            END DO 
     1186         END DO 
     1187      END DO 
     1188      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
     1189      ! 
     1190      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1191         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
     1192      END DO 
     1193!!gm  bug ? :  must be a do loop with mj0,mj1 
     1194      !  
     1195      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
     1196      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
     1197      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
     1198      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
     1199      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
     1200 
     1201      ! Control of the sign 
     1202      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
     1203      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
     1204      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
     1205      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
     1206      
     1207      ! Compute gdep3w_0 (vertical sum of e3w) 
     1208      IF ( ln_isfcav ) THEN ! if cavity 
     1209         WHERE (misfdep == 0) misfdep = 1 
     1210         DO jj = 1,jpj 
     1211            DO ji = 1,jpi 
     1212               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
     1213               DO jk = 2, misfdep(ji,jj) 
     1214                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1215               END DO 
     1216               IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
     1217               DO jk = misfdep(ji,jj) + 1, jpk 
     1218                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1219               END DO 
     1220            END DO 
     1221         END DO 
     1222      ELSE ! no cavity 
     1223         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
     1224         DO jk = 2, jpk 
     1225            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     1226         END DO 
     1227      END IF 
     1228      !                                               ! ================= ! 
     1229      IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
     1230         !                                            ! ================= ! 
     1231         DO jj = 1,jpj 
     1232            DO ji = 1, jpi 
     1233               ik = MAX( mbathy(ji,jj), 1 ) 
     1234               zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
     1235               zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
     1236               zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
     1237               zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
     1238               zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
     1239               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
     1240            END DO 
     1241         END DO 
     1242         WRITE(numout,*) 
     1243         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1244         WRITE(numout,*) 
     1245         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1246         WRITE(numout,*) 
     1247         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1248         WRITE(numout,*) 
     1249         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1250         WRITE(numout,*) 
     1251         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1252         WRITE(numout,*) 
     1253         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1254      ENDIF   
     1255      ! 
     1256      CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
     1257      ! 
     1258      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
     1259      ! 
     1260   END SUBROUTINE zgr_zps 
     1261 
     1262   SUBROUTINE zgr_isf 
     1263      !!---------------------------------------------------------------------- 
     1264      !!                    ***  ROUTINE zgr_isf  *** 
     1265      !!    
     1266      !! ** Purpose :   check the bathymetry in levels 
     1267      !!    
     1268      !! ** Method  :   THe water column have to contained at least 2 cells 
     1269      !!                Bathymetry and isfdraft are modified (dig/close) to respect 
     1270      !!                this criterion. 
     1271      !!                  
     1272      !!    
     1273      !! ** Action  : - test compatibility between isfdraft and bathy  
     1274      !!              - bathy and isfdraft are modified 
     1275      !!---------------------------------------------------------------------- 
     1276      !!    
    9671277      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    9681278      INTEGER  ::   ik, it           ! temporary integers 
    … …  
    9751285      REAL(wp) ::   zdiff            ! temporary scalar 
    9761286      REAL(wp) ::   zrefdep          ! temporary scalar 
    977       REAL(wp) ::   zbathydiff, zrisfdepdiff  
    978       REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 3D workspace (ISH) 
    979       INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep   ! 3D workspace (ISH) 
    980       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
     1287      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar 
     1288      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
     1289      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
    9811290      !!--------------------------------------------------------------------- 
    9821291      ! 
    983       IF( nn_timing == 1 )  CALL timing_start('zgr_zps') 
    984       ! 
    985       CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     1292      IF( nn_timing == 1 )  CALL timing_start('zgr_isf') 
     1293      ! 
    9861294      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 
    987       CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep) 
    988       ! 
    989       IF(lwp) WRITE(numout,*) 
    990       IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps' 
    991       IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
    992       IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
    993  
    994       ll_print = .FALSE.                   ! Local variable for debugging 
    995        
    996       IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
    997          WRITE(numout,*) 
    998          WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
    999          CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
    1000       ENDIF 
    1001  
    1002       ! bathymetry in level (from bathy_meter) 
    1003       ! =================== 
    1004       zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 
    1005       bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
    1006       WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
    1007       ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level 
    1008       END WHERE 
    1009  
    1010       ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
    1011       ! find the number of ocean levels such that the last level thickness 
    1012       ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 
    1013       ! e3t_1d is the reference level thickness 
    1014       DO jk = jpkm1, 1, -1 
    1015          zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1016          WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
    1017       END DO 
     1295      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 
     1296 
     1297 
    10181298      ! (ISF) compute misfdep 
    10191299      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
    … …  
    10591339            misfdep(jpi,:) = misfdep(  2  ,:)  
    10601340         ENDIF 
    1061   
     1341 
    10621342         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
    10631343            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    10641344            mbathy(jpi,:) = mbathy(  2  ,:) 
    10651345         ENDIF 
    1066   
     1346 
    10671347         ! split last cell if possible (only where water column is 2 cell or less) 
    10681348         DO jk = jpkm1, 1, -1 
    … …  
    10821362            END WHERE 
    10831363         END DO 
    1084   
     1364 
    10851365  
    10861366 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
    … …  
    13631643               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
    13641644               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1365                IF( ibtest == 0 ) THEN 
     1645               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 
    13661646                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
    13671647               END IF 
    … …  
    14791759      ENDIF  
    14801760 
    1481       ! Scale factors and depth at T- and W-points 
    1482       DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
    1483          gdept_0(:,:,jk) = gdept_1d(jk) 
    1484          gdepw_0(:,:,jk) = gdepw_1d(jk) 
    1485          e3t_0  (:,:,jk) = e3t_1d  (jk) 
    1486          e3w_0  (:,:,jk) = e3w_1d  (jk) 
    1487       END DO 
    1488       !  
    1489       DO jj = 1, jpj 
    1490          DO ji = 1, jpi 
    1491             ik = mbathy(ji,jj) 
    1492             IF( ik > 0 ) THEN               ! ocean point only 
    1493                ! max ocean level case 
    1494                IF( ik == jpkm1 ) THEN 
    1495                   zdepwp = bathy(ji,jj) 
    1496                   ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
    1497                   ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
    1498                   e3t_0(ji,jj,ik  ) = ze3tp 
    1499                   e3t_0(ji,jj,ik+1) = ze3tp 
    1500                   e3w_0(ji,jj,ik  ) = ze3wp 
    1501                   e3w_0(ji,jj,ik+1) = ze3tp 
    1502                   gdepw_0(ji,jj,ik+1) = zdepwp 
    1503                   gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
    1504                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
    1505                   ! 
    1506                ELSE                         ! standard case 
    1507                   IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
    1508                   ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    1509                   ENDIF 
    1510 !gm Bug?  check the gdepw_1d 
    1511                   !       ... on ik 
    1512                   gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
    1513                      &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    1514                      &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1515                   e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    1516                      &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    1517                   e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    1518                      &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    1519                   !       ... on ik+1 
    1520                   e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1521                   e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1522                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    1523                ENDIF 
    1524             ENDIF 
    1525          END DO 
    1526       END DO 
    1527       ! 
    1528       it = 0 
    1529       DO jj = 1, jpj 
    1530          DO ji = 1, jpi 
    1531             ik = mbathy(ji,jj) 
    1532             IF( ik > 0 ) THEN               ! ocean point only 
    1533                e3tp (ji,jj) = e3t_0(ji,jj,ik) 
    1534                e3wp (ji,jj) = e3w_0(ji,jj,ik) 
    1535                ! test 
    1536                zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
    1537                IF( zdiff <= 0._wp .AND. lwp ) THEN  
    1538                   it = it + 1 
    1539                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
    1540                   WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
    1541                   WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
    1542                   WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
    1543                ENDIF 
    1544             ENDIF 
    1545          END DO 
    1546       END DO 
    1547       ! 
    1548       ! (ISF) Definition of e3t, u, v, w for ISF case 
    1549       DO jj = 1, jpj  
    1550          DO ji = 1, jpi  
    1551             ik = misfdep(ji,jj)  
    1552             IF( ik > 1 ) THEN               ! ice shelf point only  
    1553                IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
    1554                gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
    1555 !gm Bug?  check the gdepw_0  
    1556                !       ... on ik  
    1557                gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
    1558                   &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
    1559                   &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
    1560                e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
    1561                e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
    1562  
    1563                IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
    1564                   e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    1565                 ENDIF  
    1566                !       ... on ik / ik-1  
    1567                e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    1568                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    1569 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
    1570                gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
    1571             ENDIF  
    1572          END DO  
    1573       END DO  
    1574       !  
    1575       it = 0  
    1576       DO jj = 1, jpj  
    1577          DO ji = 1, jpi  
    1578             ik = misfdep(ji,jj)  
    1579             IF( ik > 1 ) THEN               ! ice shelf point only  
    1580                e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
    1581                e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
    1582                ! test  
    1583                zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
    1584                IF( zdiff <= 0. .AND. lwp ) THEN   
    1585                   it = it + 1  
    1586                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
    1587                   WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
    1588                   WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
    1589                   WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
    1590                ENDIF  
    1591             ENDIF  
    1592          END DO  
    1593       END DO  
    1594       ! END (ISF) 
    1595  
    1596       ! Scale factors and depth at U-, V-, UW and VW-points 
    1597       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1598          e3u_0 (:,:,jk) = e3t_1d(jk) 
    1599          e3v_0 (:,:,jk) = e3t_1d(jk) 
    1600          e3uw_0(:,:,jk) = e3w_1d(jk) 
    1601          e3vw_0(:,:,jk) = e3w_1d(jk) 
    1602       END DO 
    1603       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
    1604          DO jj = 1, jpjm1 
    1605             DO ji = 1, fs_jpim1   ! vector opt. 
    1606                e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
    1607                e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
    1608                e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
    1609                e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
    1610             END DO 
    1611          END DO 
    1612       END DO 
    1613       ! (ISF) define e3uw 
    1614       DO jk = 2,jpk                           
    1615          DO jj = 1, jpjm1  
    1616             DO ji = 1, fs_jpim1   ! vector opt.  
    1617                e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) & 
    1618                  &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) ) 
    1619                e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) & 
    1620                  &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) ) 
    1621             END DO  
    1622          END DO  
    1623       END DO 
    1624       !End (ISF) 
    1625        
    1626       CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    1627       CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    1628       ! 
    1629       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1630          WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
    1631          WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
    1632          WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
    1633          WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
    1634       END DO 
    1635        
    1636       ! Scale factor at F-point 
    1637       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1638          e3f_0(:,:,jk) = e3t_1d(jk) 
    1639       END DO 
    1640       DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
    1641          DO jj = 1, jpjm1 
    1642             DO ji = 1, fs_jpim1   ! vector opt. 
    1643                e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
    1644             END DO 
    1645          END DO 
    1646       END DO 
    1647       CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
    1648       ! 
    1649       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1650          WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
    1651       END DO 
    1652 !!gm  bug ? :  must be a do loop with mj0,mj1 
    1653       !  
    1654       e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
    1655       e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
    1656       e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
    1657       e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
    1658       e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
    1659  
    1660       ! Control of the sign 
    1661       IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
    1662       IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
    1663       IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
    1664       IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
    1665       
    1666       ! Compute gdep3w_0 (vertical sum of e3w) 
    1667       WHERE (misfdep == 0) misfdep = 1 
    1668       DO jj = 1,jpj 
    1669          DO ji = 1,jpi 
    1670             gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
    1671             DO jk = 2, misfdep(ji,jj) 
    1672                gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1673             END DO 
    1674             IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
    1675             DO jk = misfdep(ji,jj) + 1, jpk 
    1676                gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1677             END DO 
    1678         END DO 
    1679       END DO 
    1680       !                                               ! ================= ! 
    1681       IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
    1682          !                                            ! ================= ! 
    1683          DO jj = 1,jpj 
    1684             DO ji = 1, jpi 
    1685                ik = MAX( mbathy(ji,jj), 1 ) 
    1686                zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
    1687                zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
    1688                zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
    1689                zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
    1690                zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
    1691                zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
    1692             END DO 
    1693          END DO 
    1694          WRITE(numout,*) 
    1695          WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1696          WRITE(numout,*) 
    1697          WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1698          WRITE(numout,*) 
    1699          WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1700          WRITE(numout,*) 
    1701          WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1702          WRITE(numout,*) 
    1703          WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1704          WRITE(numout,*) 
    1705          WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1706       ENDIF   
    1707       ! 
    1708       CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
    17091761      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 
    17101762      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 
    1711       ! 
    1712       IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
    1713       ! 
    1714    END SUBROUTINE zgr_zps 
     1763 
     1764      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf') 
     1765       
     1766   END SUBROUTINE 
    17151767 
    17161768   SUBROUTINE zgr_sco 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r4990 r5120  
    137137         CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )        ! before potential and in situ densities 
    138138#if ! defined key_c1d 
    139          IF( ln_zps )    CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    140             &                                      rhd, gru , grv, aru, arv, gzu, gzv, ge3ru, ge3rv,  &             ! 
    141             &                                      gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     139         IF( ln_zps .AND. .NOT. ln_isfcav)                                 & 
     140            &            CALL zps_hde    ( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
     141            &                                            rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     142         IF( ln_zps .AND.       ln_isfcav)                                 & 
     143            &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     144            &                                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     145            &                                     gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    142146#endif 
    143147         !    

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

    r4990 r5120  
    1717   !!            3.3  ! 2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
    1818   !!             -   ! 2010-10  (R. Furner, G. Madec) runoff and cla added directly here 
     19   !!            3.6  ! 2014-11  (P. Mathiot)          isf            added directly here 
    1920   !!---------------------------------------------------------------------- 
    2021 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90

    r4990 r5120  
    127127      IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec )   & 
    128128          CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' ) 
     129      IF( ln_dynzad_zts .AND. ln_isfcav )   & 
     130          CALL ctl_stop( 'Sub timestepping of vertical advection does not work with ln_isfcav = .TRUE.' ) 
    129131 
    130132      !                               ! Set nadv 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90

    r4990 r5120  
    8080              ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) 
    8181              va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) 
    82                
    83               ! (ISF) stability criteria for top friction 
    84               ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels 
    85               ikbv = mikv(ji,jj) 
    86               ! 
    87               ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
    88               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  tfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) & 
    89                  &             * (1.-umask(ji,jj,1)) 
    90               va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  tfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) & 
    91                  &             * (1.-vmask(ji,jj,1)) 
    92               ! (ISF) 
    93                
    9482           END DO 
    9583        END DO 
     84         
     85        IF ( ln_isfcav ) THEN 
     86           DO jj = 2, jpjm1 
     87              DO ji = 2, jpim1 
     88                 ! (ISF) stability criteria for top friction 
     89                 ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels 
     90                 ikbv = mikv(ji,jj) 
     91                 ! 
     92                 ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
     93                 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  tfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) & 
     94                    &             * (1.-umask(ji,jj,1)) 
     95                 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  tfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) & 
     96                    &             * (1.-vmask(ji,jj,1)) 
     97                 ! (ISF) 
     98              END DO 
     99           END DO 
     100        END IF 
    96101 
    97102        ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r4990 r5120  
    1616   !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates 
    1717   !!                 !           (A. Coward) suppression of hel, wdj and rot options 
     18   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity 
    1819   !!---------------------------------------------------------------------- 
    1920 
    … …  
    2526   !!       hpg_zps  : z-coordinate plus partial steps (interpolation) 
    2627   !!       hpg_sco  : s-coordinate (standard jacobian formulation) 
     28   !!       hpg_isf  : s-coordinate (sco formulation) adapted to ice shelf 
    2729   !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial) 
    2830   !!       hpg_prj  : s-coordinate (Pressure Jacobian with Cubic polynomial) 
    … …  
    5557   LOGICAL , PUBLIC ::   ln_hpg_djc      !: s-coordinate (Density Jacobian with Cubic polynomial) 
    5658   LOGICAL , PUBLIC ::   ln_hpg_prj      !: s-coordinate (Pressure Jacobian scheme) 
     59   LOGICAL , PUBLIC ::   ln_hpg_isf      !: s-coordinate similar to sco modify for isf 
    5760   LOGICAL , PUBLIC ::   ln_dynhpg_imp   !: semi-implicite hpg flag 
    5861 
    … …  
    97100      CASE (  3 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
    98101      CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme) 
     102      CASE (  5 )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf 
    99103      END SELECT 
    100104      ! 
    … …  
    128132      !! 
    129133      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     & 
    130          &                 ln_hpg_djc, ln_hpg_prj, ln_dynhpg_imp 
     134         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, ln_dynhpg_imp 
    131135      !!---------------------------------------------------------------------- 
    132136      ! 
    … …  
    148152         WRITE(numout,*) '      z-coord. - partial steps (interpolation)          ln_hpg_zps    = ', ln_hpg_zps 
    149153         WRITE(numout,*) '      s-coord. (standard jacobian formulation)          ln_hpg_sco    = ', ln_hpg_sco 
     154         WRITE(numout,*) '      s-coord. (standard jacobian formulation) for isf  ln_hpg_isf    = ', ln_hpg_isf 
    150155         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc 
    151156         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj 
    … …  
    158163                           & either  ln_hpg_sco or  ln_hpg_prj instead') 
    159164      ! 
    160       IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) )   & 
     165      IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )   & 
    161166         &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 
    162167                           & the standard jacobian formulation hpg_sco or & 
    163168                           & the pressure jacobian formulation hpg_prj') 
     169 
     170      IF(       ln_hpg_isf .AND. .NOT. ln_isfcav )   & 
     171         &   CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' ) 
     172      IF( .NOT. ln_hpg_isf .AND.       ln_isfcav )   & 
     173         &   CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' ) 
    164174      ! 
    165175      !                               ! Set nhpg from ln_hpg_... flags 
    … …  
    169179      IF( ln_hpg_djc )   nhpg = 3 
    170180      IF( ln_hpg_prj )   nhpg = 4 
     181      IF( ln_hpg_isf )   nhpg = 5 
    171182      ! 
    172183      !                               ! Consistency check 
    … …  
    177188      IF( ln_hpg_djc )   ioptio = ioptio + 1 
    178189      IF( ln_hpg_prj )   ioptio = ioptio + 1 
     190      IF( ln_hpg_isf )   ioptio = ioptio + 1 
    179191      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 
    180       IF( (ln_hpg_zco .OR. ln_hpg_zps .OR. ln_hpg_djc .OR. ln_hpg_prj ) .AND. nn_isf .NE. 0 )   & 
    181           &  CALL ctl_stop( 'Only hpg_sco has been corrected to work with ice shelf cavity.' ) 
     192      !  
     193      ! initialisation of ice load 
     194      riceload(:,:)=0.0 
    182195      ! 
    183196   END SUBROUTINE dyn_hpg_init 
    … …  
    345358   END SUBROUTINE hpg_zps 
    346359 
    347  
    348360   SUBROUTINE hpg_sco( kt ) 
    349361      !!--------------------------------------------------------------------- 
    … …  
    366378      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    367379      !! 
     380      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
     381      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
     382      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
     383      !!---------------------------------------------------------------------- 
     384      ! 
     385      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     386      ! 
     387      IF( kt == nit000 ) THEN 
     388         IF(lwp) WRITE(numout,*) 
     389         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 
     390         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used' 
     391      ENDIF 
     392 
     393      ! Local constant initialization 
     394      zcoef0 = - grav * 0.5_wp 
     395      ! To use density and not density anomaly 
     396      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume 
     397      ELSE                 ;     znad = 0._wp         ! Fixed volume 
     398      ENDIF 
     399 
     400      ! Surface value 
     401      DO jj = 2, jpjm1 
     402         DO ji = fs_2, fs_jpim1   ! vector opt. 
     403            ! hydrostatic pressure gradient along s-surfaces 
     404            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   & 
     405               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
     406            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   & 
     407               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
     408            ! s-coordinate pressure gradient correction 
     409            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
     410               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
     411            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
     412               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     413            ! add to the general momentum trend 
     414            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     415            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
     416         END DO 
     417      END DO 
     418 
     419      ! interior value (2=<jk=<jpkm1) 
     420      DO jk = 2, jpkm1 
     421         DO jj = 2, jpjm1 
     422            DO ji = fs_2, fs_jpim1   ! vector opt. 
     423               ! hydrostatic pressure gradient along s-surfaces 
     424               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
     425                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
     426                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
     427               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
     428                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
     429                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
     430               ! s-coordinate pressure gradient correction 
     431               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     432                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
     433               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     434                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     435               ! add to the general momentum trend 
     436               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     437               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     438            END DO 
     439         END DO 
     440      END DO 
     441      ! 
     442      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     443      ! 
     444   END SUBROUTINE hpg_sco 
     445 
     446   SUBROUTINE hpg_isf( kt ) 
     447      !!--------------------------------------------------------------------- 
     448      !!                  ***  ROUTINE hpg_sco  *** 
     449      !! 
     450      !! ** Method  :   s-coordinate case. Jacobian scheme. 
     451      !!      The now hydrostatic pressure gradient at a given level, jk, 
     452      !!      is computed by taking the vertical integral of the in-situ 
     453      !!      density gradient along the model level from the suface to that 
     454      !!      level. s-coordinates (ln_sco): a corrective term is added 
     455      !!      to the horizontal pressure gradient : 
     456      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
     457      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
     458      !!      add it to the general momentum trend (ua,va). 
     459      !!         ua = ua - 1/e1u * zhpi 
     460      !!         va = va - 1/e2v * zhpj 
     461      !!      iceload is added and partial cell case are added to the top and bottom 
     462      !!       
     463      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     464      !!---------------------------------------------------------------------- 
     465      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     466      !! 
    368467      INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j                 ! dummy loop indices 
    369468      REAL(wp) ::   zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1             ! temporary scalars 
    … …  
    379478     IF( kt == nit000 ) THEN 
    380479         IF(lwp) WRITE(numout,*) 
    381          IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 
     480         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 
    382481         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used' 
    383482      ENDIF 
    … …  
    565664!================================================================================== 
    566665 
    567 # if defined key_vectopt_loop 
    568          jj = 1 
    569          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    570 # else 
    571666      DO jj = 2, jpjm1 
    572667         DO ji = 2, jpim1 
    573 # endif 
    574668            iku = mbku(ji,jj) 
    575669            ikv = mbkv(ji,jj) 
    … …  
    598692               va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 
    599693            END IF 
    600 # if ! defined key_vectopt_loop 
    601          END DO 
    602 # endif 
     694         END DO 
    603695      END DO 
    604696      
    … …  
    610702      CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 
    611703      ! 
    612    END SUBROUTINE hpg_sco 
     704   END SUBROUTINE hpg_isf 
    613705 
    614706 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r4990 r5120  
    250250      IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   & 
    251251           &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 
    252       IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. nn_isf .NE. 0 )   & 
     252      IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav )   & 
    253253           &   CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 
    254254      ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r5032 r5120  
    2222   USE dom_oce         ! ocean space and time domain 
    2323   USE sbc_oce         ! surface boundary condition: ocean 
     24   USE sbcisf          ! ice shelf variable (fwfisf) 
    2425   USE dynspg_oce      ! surface pressure gradient variables 
    2526   USE phycst          ! physical constants 
    … …  
    453454      !                                         ! Surface net water flux and rivers 
    454455      IF (ln_bt_fw) THEN 
    455          zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) ) 
     456         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) ) 
    456457      ELSE 
    457          zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)) 
     458         zssh_frc(:,:) = zraur * z1_2 * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   & 
     459                &                        + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) )       ) 
    458460      ENDIF 
    459461#if defined key_asminc 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r4990 r5120  
    9595         END DO    
    9696      END DO 
    97       DO jj = 2, jpjm1              ! Surface and bottom values set to zero 
    98          DO ji = fs_2, fs_jpim1           ! vector opt. 
    99             zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp 
    100             zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp 
    101             zwuw(ji,jj,jpk) = 0._wp 
    102             zwvw(ji,jj,jpk) = 0._wp 
    103          END DO   
    104       END DO 
     97      ! 
     98      ! Surface and bottom advective fluxes set to zero 
     99      IF ( ln_isfcav ) THEN 
     100         DO jj = 2, jpjm1 
     101            DO ji = fs_2, fs_jpim1           ! vector opt. 
     102               zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp 
     103               zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp 
     104               zwuw(ji,jj,jpk) = 0._wp 
     105               zwvw(ji,jj,jpk) = 0._wp 
     106            END DO 
     107         END DO 
     108      ELSE 
     109         DO jj = 2, jpjm1         
     110            DO ji = fs_2, fs_jpim1           ! vector opt. 
     111               zwuw(ji,jj, 1 ) = 0._wp 
     112               zwvw(ji,jj, 1 ) = 0._wp 
     113               zwuw(ji,jj,jpk) = 0._wp 
     114               zwvw(ji,jj,jpk) = 0._wp 
     115            END DO   
     116         END DO 
     117      END IF 
    105118 
    106119      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points 
    … …  
    196209         END DO 
    197210      END DO 
    198  
    199       DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero 
     211      ! 
     212      ! Surface and bottom advective fluxes set to zero 
     213      DO jj = 2, jpjm1         
    200214         DO ji = fs_2, fs_jpim1           ! vector opt. 
    201             zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp 
    202             zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp 
     215            zwuw(ji,jj, 1 ) = 0._wp 
     216            zwvw(ji,jj, 1 ) = 0._wp 
    203217            zwuw(ji,jj,jpk) = 0._wp 
    204218            zwvw(ji,jj,jpk) = 0._wp 
    … …  
    228242            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point 
    229243               DO ji = fs_2, fs_jpim1        ! vector opt. 
    230                   zwuw(ji,jj,jk) = ( zww(ji+1,jj  ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) 
    231                   zwvw(ji,jj,jk) = ( zww(ji  ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) 
     244                  zwuw(ji,jj,jk) = ( zww(ji+1,jj  ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) !* wumask(ji,jj,jk) 
     245                  zwvw(ji,jj,jk) = ( zww(ji  ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) !* wvmask(ji,jj,jk) 
    232246               END DO   
    233247            END DO    

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r4990 r5120  
    105105               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
    106106               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
    107                ikbu = miku(ji,jj)       ! ocean top level at u- and v-points  
    108                ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
    109                IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu) 
    110                IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv) 
    111             END DO 
    112          END DO 
     107            END DO 
     108         END DO 
     109         IF ( ln_isfcav ) THEN 
     110            DO jj = 2, jpjm1 
     111               DO ji = 2, jpim1 
     112                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points  
     113                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
     114                  IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu) 
     115                  IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv) 
     116               END DO 
     117            END DO 
     118         END IF 
    113119      ENDIF 
    114120 
    … …  
    145151               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
    146152               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 
    147                ikbu = miku(ji,jj)         ! top ocean level at u- and v-points  
    148                ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
    149                ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu) 
    150                ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv) 
    151                ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
    152                va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va 
    153             END DO 
    154          END DO 
     153            END DO 
     154         END DO 
     155         IF ( ln_isfcav ) THEN 
     156            DO jj = 2, jpjm1         
     157               DO ji = fs_2, fs_jpim1   ! vector opt. 
     158                  ikbu = miku(ji,jj)         ! top ocean level at u- and v-points  
     159                  ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
     160                  ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu) 
     161                  ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv) 
     162                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
     163                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va 
     164               END DO 
     165            END DO 
     166         END IF 
    155167      ENDIF 
    156168#endif 
    … …  
    167179               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point 
    168180               zcoef = - p2dt / ze3ua       
    169                zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
    170                zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk) 
    171                zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)  
    172                zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1) 
    173                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
     181               zzwi          = zcoef * avmu  (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
     182               zwi(ji,jj,jk) = zzwi  * wumask(ji,jj,jk  ) 
     183               zzws          = zcoef * avmu  (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)  
     184               zws(ji,jj,jk) = zzws  * wumask(ji,jj,jk+1) 
     185               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    174186            END DO 
    175187         END DO 
    … …  
    198210      ! 
    199211      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    200       DO jj = 2, jpjm1    
    201          DO ji = fs_2, fs_jpim1   ! vector opt. 
    202             DO jk = miku(ji,jj)+1, jpkm1 
     212      DO jk = 2, jpkm1 
     213         DO jj = 2, jpjm1    
     214            DO ji = fs_2, fs_jpim1   ! vector opt. 
    203215               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    204216            END DO 
    … …  
    208220      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    209221         DO ji = fs_2, fs_jpim1   ! vector opt. 
    210             ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,miku(ji,jj)) + r_vvl   * fse3u_a(ji,jj,miku(ji,jj))  
    211222#if defined key_dynspg_ts 
    212             ua(ji,jj,miku(ji,jj)) = ua(ji,jj,miku(ji,jj)) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    213                &                                      / ( ze3ua * rau0 )  
     223            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1)  
     224            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     225               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
    214226#else 
    215             ua(ji,jj,miku(ji,jj)) = ub(ji,jj,miku(ji,jj)) & 
    216                &                   + p2dt *(ua(ji,jj,miku(ji,jj)) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    217                &                                  / ( fse3u(ji,jj,miku(ji,jj)) * rau0     ) )  
    218 #endif 
    219             DO jk = miku(ji,jj)+1, jpkm1 
     227            ua(ji,jj,1) = ub(ji,jj,1) & 
     228               &                   + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     229               &                                      / ( fse3u(ji,jj,1) * rau0     ) * umask(ji,jj,1) )  
     230#endif 
     231         END DO 
     232      END DO 
     233      DO jk = 2, jpkm1 
     234         DO jj = 2, jpjm1 
     235            DO ji = fs_2, fs_jpim1 
    220236#if defined key_dynspg_ts 
    221237               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side 
    … …  
    231247         DO ji = fs_2, fs_jpim1   ! vector opt. 
    232248            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    233             DO jk = jpk-2, miku(ji,jj), -1 
     249         END DO 
     250      END DO 
     251      DO jk = jpk-2, 1, -1 
     252         DO jj = 2, jpjm1 
     253            DO ji = fs_2, fs_jpim1 
    234254               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    235255            END DO 
    … …  
    260280               zcoef = - p2dt / ze3va 
    261281               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
    262                zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk) 
     282               zwi(ji,jj,jk) =  zzwi * wvmask(ji,jj,jk) 
    263283               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
    264                zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1) 
    265                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
     284               zws(ji,jj,jk) =  zzws * wvmask(ji,jj,jk+1) 
     285               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    266286            END DO 
    267287         END DO 
    … …  
    290310      ! 
    291311      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    292       DO jj = 2, jpjm1    
    293          DO ji = fs_2, fs_jpim1   ! vector opt. 
    294             DO jk = mikv(ji,jj)+1, jpkm1         
     312      DO jk = 2, jpkm1         
     313         DO jj = 2, jpjm1    
     314            DO ji = fs_2, fs_jpim1   ! vector opt. 
    295315               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    296316            END DO 
    … …  
    300320      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    301321         DO ji = fs_2, fs_jpim1   ! vector opt. 
    302             ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,mikv(ji,jj)) + r_vvl   * fse3v_a(ji,jj,mikv(ji,jj))  
    303322#if defined key_dynspg_ts             
    304             va(ji,jj,mikv(ji,jj)) = va(ji,jj,mikv(ji,jj)) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     323            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1)  
     324            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    305325               &                                      / ( ze3va * rau0 )  
    306326#else 
    307             va(ji,jj,mikv(ji,jj)) = vb(ji,jj,mikv(ji,jj)) & 
    308                &                   + p2dt *(va(ji,jj,mikv(ji,jj)) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    309                &                                                       / ( fse3v(ji,jj,mikv(ji,jj)) * rau0     )  ) 
    310 #endif 
    311             DO jk = mikv(ji,jj)+1, jpkm1 
     327            va(ji,jj,1) = vb(ji,jj,1) & 
     328               &                   + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     329               &                                                       / ( fse3v(ji,jj,1) * rau0     )  ) 
     330#endif 
     331         END DO 
     332      END DO 
     333      DO jk = 2, jpkm1 
     334         DO jj = 2, jpjm1 
     335            DO ji = fs_2, fs_jpim1   ! vector opt. 
    312336#if defined key_dynspg_ts 
    313337               zrhs = va(ji,jj,jk)   ! zrhs=right hand side 
    … …  
    323347         DO ji = fs_2, fs_jpim1   ! vector opt. 
    324348            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    325             DO jk = jpk-2, mikv(ji,jj), -1 
     349         END DO 
     350      END DO 
     351      DO jk = jpk-2, 1, -1 
     352         DO jj = 2, jpjm1 
     353            DO ji = fs_2, fs_jpim1 
    326354               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    327355            END DO 
    … …  
    349377              avmu(ji,jj,ikbu+1) = 0.e0 
    350378              avmv(ji,jj,ikbv+1) = 0.e0 
    351               ikbu = miku(ji,jj)         ! ocean top level at u- and v-points  
    352               ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
    353               IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0 
    354               IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0 
    355379           END DO 
    356380        END DO 
     381        IF (ln_isfcav) THEN 
     382           DO jj = 2, jpjm1 
     383              DO ji = 2, jpim1 
     384                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points  
     385                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
     386                 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0 
     387                 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0 
     388              END DO 
     389           END DO 
     390        END IF 
    357391      ENDIF 
    358392      ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r5016 r5120  
    142142            DO jj = 1, jpjm1 
    143143               DO ji = 1, jpim1 
    144 ! IF should be useless check zpshde (PM) 
    145                IF ( mbku(ji,jj) > 1 ) zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
    146                IF ( mbkv(ji,jj) > 1 ) zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 
     144                  zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
     145                  zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 
     146               END DO 
     147            END DO 
     148         ENDIF 
     149         IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level 
     150            DO jj = 1, jpjm1 
     151               DO ji = 1, jpim1 
    147152               IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj)  
    148153               IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 
    … …  
    151156         ENDIF 
    152157         ! 
    153          zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2) 
    154          DO jk = 1, jpkm1 
     158         !==   Local vertical density gradient at T-point   == !   (evaluated from N^2) 
     159         ! interior value 
     160         DO jk = 2, jpkm1 
    155161            !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 
    156162            !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0 
    … …  
    162168         END DO 
    163169         ! surface initialisation  
    164          DO jj = 1, jpjm1 
    165             DO ji = 1, jpim1 
    166               zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp 
    167             END DO 
    168          END DO 
     170         zdzr(:,:,1) = 0._wp  
     171         IF ( ln_isfcav ) THEN 
     172            ! if isf need to overwrite the interior value at at the first ocean point 
     173            DO jj = 1, jpjm1 
     174               DO ji = 1, jpim1 
     175                  zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp 
     176               END DO 
     177            END DO 
     178         END IF 
    169179         ! 
    170180         !                          !==   Slopes just below the mixed layer   ==! 
    … …  
    175185         ! ===========================      | vslp = d/dj( prd ) / d/dz( prd ) 
    176186         ! 
    177          DO jj = 2, jpjm1 
    178             DO ji = fs_2, fs_jpim1   ! vector opt. 
    179                IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = hmlpt(ji  ,jj) 
    180                IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = hmlpt(ji+1,jj) 
    181                IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj), hmlpt(ji+1,jj)) 
    182                IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = hmlpt(ji  ,jj) 
    183                IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = hmlpt(ji,jj+1) 
    184                IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji,jj+1)) 
     187         IF ( ln_isfcav ) THEN 
     188            DO jj = 2, jpjm1 
     189               DO ji = fs_2, fs_jpim1   ! vector opt. 
     190                  IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp) 
     191                  IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji+1,jj  ),                   5._wp) 
     192                  IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji+1,jj  ), 5._wp) 
     193                  IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp) 
     194                  IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj+1),                   5._wp) 
     195                  IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji  ,jj+1), 5._wp) 
     196               ENDDO 
    185197            ENDDO 
    186          ENDDO 
     198         ELSE 
     199            DO jj = 2, jpjm1 
     200               DO ji = fs_2, fs_jpim1   ! vector opt. 
     201                  zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp) 
     202                  zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp) 
     203               ENDDO 
     204            ENDDO 
     205         END IF 
    187206         DO jk = 2, jpkm1                            !* Slopes at u and v points 
    188207            DO jj = 2, jpjm1 
    … …  
    198217                  zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  ) 
    199218                  !                                      ! uslp and vslp output in zwz and zww, resp. 
    200                   zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )  
    201                   zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )  
     219                  zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj  ,jk) ) 
     220                  zfj = MAX( omlmask(ji,jj,jk), omlmask(ji  ,jj+1,jk) ) 
    202221                  ! thickness of water column between surface and level k at u/v point 
    203                   zdepu = 0.5_wp * (( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                   & 
    204                              - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) )  & 
    205                              - fse3u(ji,jj,miku(ji,jj))                                         ) 
    206                   zdepv = 0.5_wp * (( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                   & 
    207                              - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) & 
    208                              - fse3v(ji,jj,mikv(ji,jj))                                         ) 
    209                   zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                              & 
    210                      &                 + zfi  * uslpml(ji,jj)                                                     & 
    211                      &                        * zdepu / MAX( zhmlpu(ji,jj), 5._wp ) 
    212                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * umask(ji,jj,jk) * umask(ji,jj,jk-1) 
    213                   zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                              & 
    214                      &                 + zfj  * vslpml(ji,jj)                                                     & 
    215                      &                        * zdepv / MAX( zhmlpv(ji,jj), 5._wp )  
    216                   zww(ji,jj,jk) = zww(ji,jj,jk) * vmask(ji,jj,jk) * vmask(ji,jj,jk-1) 
     222                  zdepu = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                              & 
     223                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) ) - fse3u(ji,jj,miku(ji,jj)) ) 
     224                  zdepv = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                              & 
     225                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji  ,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) ) 
     226                  ! 
     227                  zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                          & 
     228                     &                 + zfi  * uslpml(ji,jj) * zdepu / zhmlpu(ji,jj) 
     229                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk) 
     230                  zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                          & 
     231                     &                 + zfj  * vslpml(ji,jj) * zdepv / zhmlpv(ji,jj)  
     232                  zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk) 
    217233                   
    218234                  
    … …  
    266282                  uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   & 
    267283                     &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp   & 
    268                      &                            *   umask(ji,jj,jk-1) !* umask(ji,jj,jk) * umask(ji,jj,jk+1) 
     284                     &                            *   umask(ji,jj,jk-1) 
    269285                  vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   & 
    270286                     &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp   & 
    271                      &                            *   vmask(ji,jj,jk-1) !* vmask(ji,jj,jk) * vmask(ji,jj,jk+1) 
     287                     &                            *   vmask(ji,jj,jk-1) 
    272288               END DO 
    273289            END DO 
    … …  
    282298               DO ji = fs_2, fs_jpim1   ! vector opt. 
    283299                  !                                  !* Local vertical density gradient evaluated from N^2 
    284                   zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     300                  zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * wmask(ji,jj,jk) 
    285301                  !                                  !* Slopes at w point 
    286302                  !                                        ! i- & j-gradient of density at w-points 
    … …  
    298314                  zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  ) 
    299315                  !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 
    300                   zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )    ! zfk=1 in the ML otherwise zfk=0 
     316                  zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0 
    301317                  zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp ) 
    302318                  zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) & 
    303                      &            + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     319                     &            + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk) 
    304320                  zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) & 
    305                      &            + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     321                     &            + zck * wslpjml(ji,jj) * zfk  ) * wmask(ji,jj,jk) 
    306322 
    307323!!gm  modif to suppress omlmask....  (as in Griffies operator) 
    … …  
    356372                  zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   & 
    357373                     &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 
    358                   wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * tmask(ji,jj,jk-1) * tmask(ji,jj,jk) 
    359                   wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * tmask(ji,jj,jk-1) * tmask(ji,jj,jk) 
     374                  wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * wmask(ji,jj,jk) 
     375                  wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * wmask(ji,jj,jk) 
    360376               END DO 
    361377            END DO 
    … …  
    423439                  vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)  
    424440                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) & 
    425                     &                              * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * 0.5  
     441                    &                              * wmask(ji,jj,jk) * 0.5  
    426442                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) & 
    427                     &                              * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * 0.5  
     443                    &                              * wmask(ji,jj,jk) * 0.5  
    428444               END DO  
    429445            END DO  
    … …  
    736752            DO ji = 1, jpi 
    737753               ik = nmln(ji,jj) - 1 
    738                IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN   ;   omlmask(ji,jj,jk) = 1._wp 
    739                ELSE                  ;   omlmask(ji,jj,jk) = 0._wp 
     754               IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN 
     755                  omlmask(ji,jj,jk) = 1._wp 
     756               ELSE 
     757                  omlmask(ji,jj,jk) = 0._wp 
    740758               ENDIF 
    741759            END DO 
    … …  
    794812            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj )  ) 
    795813            !                        !- i- & j-slope at w-points (wslpiml, wslpjml) 
    796             wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 
    797             wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 
     814            wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik) 
     815            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik) 
    798816         END DO 
    799817      END DO 

  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r4990 r5120  
    9898   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp_tot           !: total E-P over ocean and ice                 [Kg/m2/s] 
    9999   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmflx            !: freshwater budget: freezing/melting          [Kg/m2/s] 
    100    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff   [Kg/m2/s]   
     100   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff        [Kg/m2/s]   
     101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwfisf , fwfisf_b !: ice shelf melting   [Kg/m2/s]   
    101102   !! 
    102103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  sbc_tsc, sbc_tsc_b  !: sbc content trend                      [K.m/s] jpi,jpj,jpts 
    … …  
    147148         &      sfx    (jpi,jpj) , sfx_b(jpi,jpj) , emp_tot(jpi,jpj), fmmflx(jpi,jpj), STAT=ierr(2) ) 
    148149         ! 
    149       ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     & 
    150          &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) ) 
     150      ALLOCATE( fwfisf  (jpi,jpj), rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     & 
     151         &      fwfisf_b(jpi,jpj), rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) ) 
    151152         ! 
    152153      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     & 

  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90

    r4990 r5120  
    88   !!            3.0  ! 2006-08  (G. Madec)  Surface module 
    99   !!            3.2  ! 2009-07  (C. Talandier) emp mean s spread over erp area  
     10   !!            3.6  ! 2014-11  (P. Mathiot  ) add ice shelf melting 
    1011   !!---------------------------------------------------------------------- 
    1112 
    … …  
    8889         ! 
    8990         IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' ) 
    90          ! 
    91          area = glob_sum( e1e2t(:,:) )           ! interior global domain surface 
     91         IF( kn_fwb == 3 .AND. ln_isfcav    )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' ) 
     92         ! 
     93         area = glob_sum( e1e2t(:,:) * tmask(:,:,1))           ! interior global domain surface 
     94         ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes 
     95         ! and in case of no melt, it can generate HSSW. 
    9296         ! 
    9397#if ! defined key_lim2 &&  ! defined key_lim3 && ! defined key_cice 
    … …  
    106110            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain 
    107111            zcoef = z_fwf * rcp 
    108             emp(:,:) = emp(:,:) - z_fwf  
    109             qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) ! account for change to the heat budget due to fw correction 
     112            emp(:,:) = emp(:,:) - z_fwf              * tmask(:,:,1) 
     113            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction 
    110114         ENDIF 
    111115         ! 
    … …  
    138142         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes 
    139143            zcoef = fwfold * rcp 
    140             emp(:,:) = emp(:,:) + fwfold 
    141             qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) ! account for change to the heat budget due to fw correction 
     144            emp(:,:) = emp(:,:) + fwfold             * tmask(:,:,1) 
     145            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction 
    142146         ENDIF 
    143147         ! 
    … …  
    158162            zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) ) 
    159163            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges)  
    160             z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - snwice_fmass(:,:) ) ) / area 
     164            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) - snwice_fmass(:,:) ) ) / area 
    161165            !             
    162166            IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation 

  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r4990 r5120  
    77   !! History :  3.2   !  2011-02  (C.Harris  ) Original code isf cav 
    88   !!            X.X   !  2006-02  (C. Wang   ) Original code bg03 
    9    !!            3.4   !  2013-03  (P. Mathiot) Merging 
     9   !!            3.4   !  2013-03  (P. Mathiot) Merging + parametrization 
    1010   !!---------------------------------------------------------------------- 
    1111 
    … …  
    3737 
    3838   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc    
    39    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fwfisf_b, fwfisf  !: evaporation damping   [kg/m2/s] 
    40    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf            !: net heat flux from ice shelf 
     39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf              !: net heat flux from ice shelf 
    4140   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m] 
    4241   LOGICAL , PUBLIC ::   ln_divisf                   !: flag to correct divergence  
    … …  
    309308      sbc_isf_alloc = 0       ! set to zero if no array to be allocated 
    310309      IF( .NOT. ALLOCATED( qisf ) ) THEN 
    311          ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts)              , & 
    312                &    qisf(jpi,jpj)     , fwfisf(jpi,jpj)     , fwfisf_b(jpi,jpj)   , & 
    313                &    rhisf_tbl(jpi,jpj), r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , & 
    314                &    ttbl(jpi,jpj)     , stbl(jpi,jpj)       , utbl(jpi,jpj)       , & 
    315                &    vtbl(jpi, jpj)    , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), & 
    316                &    ralpha(jpi,jpj)   , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , & 
     310         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj)   , & 
     311               &    rhisf_tbl(jpi,jpj)    , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , & 
     312               &    ttbl(jpi,jpj)         , stbl(jpi,jpj)       , utbl(jpi,jpj)       , & 
     313               &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), & 
     314               &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , & 
    317315               &    STAT= sbc_isf_alloc ) 
    318316         ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r4990 r5120  
    1313   !!            3.4  ! 2011-11  (C. Harris) CICE added as an option 
    1414   !!            3.5  ! 2012-11  (A. Coward, G. Madec) Rethink of heat, mass and salt surface fluxes 
     15   !!            3.6  ! 2014-11  (P. Mathiot, C. Harris) add ice shelves melting                     
    1516   !!---------------------------------------------------------------------- 
    1617 
    … …  
    179180         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' ) 
    180181         fwfisf  (:,:) = 0.0_wp 
     182         fwfisf_b(:,:) = 0.0_wp 
    181183      END IF 
    182184      IF( nn_ice == 0  )   fr_i(:,:) = 0.e0        ! no ice in the domain, ice fraction is always zero 

  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90

    r4990 r5120  
    6161      !!--------------------------------------------------------------------- 
    6262       
    63       !                                        !* first wet T-, U-, V- ocean level (ISF) variables (T, S, depth, velocity) 
     63      !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity) 
    6464      DO jj = 1, jpj 
    6565         DO ji = 1, jpi 
    66             zub(ji,jj)        = ub (ji,jj,miku(ji,jj)) 
    67             zvb(ji,jj)        = vb (ji,jj,mikv(ji,jj)) 
    6866            zts(ji,jj,jp_tem) = tsn(ji,jj,mikt(ji,jj),jp_tem) 
    6967            zts(ji,jj,jp_sal) = tsn(ji,jj,mikt(ji,jj),jp_sal) 
    7068         END DO 
    7169      END DO 
     70      zub(:,:)        = ub (:,:,1       ) 
     71      zvb(:,:)        = vb (:,:,1       ) 
    7272      ! 
    7373      IF( lk_vvl ) THEN 
    74          DO jj = 1, jpj 
    75             DO ji = 1, jpi 
    76                zdep(ji,jj) = fse3t_n(ji,jj,mikt(ji,jj)) 
    77             END DO 
    78          END DO 
     74         zdep(:,:) = fse3t_n(:,:,1) 
    7975      ENDIF 
    8076      !                                                   ! ---------------------------------------- ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r4990 r5120  
    206206      IF( lk_esopa         )   ioptio =          1 
    207207 
    208       IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck ) .AND. nn_isf .NE. 0 )  & 
     208      IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck .OR. ln_traadv_tvd_zts ) .AND. ln_isfcav )  & 
    209209      &   CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity') 
    210210 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r4990 r5120  
    106106      ENDIF 
    107107      ! 
    108       zwi(:,:,:) = 0.e0 ; zwz(:,:,:) = 0.e0 
     108      zwi(:,:,:) = 0.e0 ;  
    109109      ! 
    110110      !                                                          ! =========== 
    111111      DO jn = 1, kjpt                                            ! tracer loop 
    112112         !                                                       ! =========== 
    113          ! 1. Bottom value : flux set to zero 
     113         ! 1. Bottom and k=1 value : flux set to zero 
    114114         ! ---------------------------------- 
    115115         zwx(:,:,jpk) = 0.e0    ;    zwz(:,:,jpk) = 0.e0 
    116116         zwy(:,:,jpk) = 0.e0    ;    zwi(:,:,jpk) = 0.e0 
    117  
     117           
     118         zwz(:,:,1  ) = 0._wp 
    118119         ! 2. upstream advection with initial mass fluxes & intermediate update 
    119120         ! -------------------------------------------------------------------- 
    … …  
    134135 
    135136         ! upstream tracer flux in the k direction 
     137         ! Interior value 
     138         DO jk = 2, jpkm1 
     139            DO jj = 1, jpj 
     140               DO ji = 1, jpi 
     141                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
     142                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     143                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
     144               END DO 
     145            END DO 
     146         END DO 
    136147         ! Surface value 
    137148         IF( lk_vvl ) THEN    
    138             DO jj = 1, jpj 
    139                DO ji = 1, jpi 
    140                   zwz(ji,jj, mikt(ji,jj) ) = 0.e0                         ! volume variable 
    141                END DO 
    142             END DO 
     149            IF ( ln_isfcav ) THEN 
     150               DO jj = 1, jpj 
     151                  DO ji = 1, jpi 
     152                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable 
     153                  END DO 
     154               END DO 
     155            ELSE 
     156               zwz(:,:,1) = 0.e0          ! volume variable 
     157            END IF 
    143158         ELSE                 
    144             DO jj = 1, jpj 
    145                DO ji = 1, jpi 
    146                   zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    147                END DO 
    148             END DO    
     159            IF ( ln_isfcav ) THEN 
     160               DO jj = 1, jpj 
     161                  DO ji = 1, jpi 
     162                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     163                  END DO 
     164               END DO    
     165            ELSE 
     166               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface 
     167            END IF 
    149168         ENDIF 
    150          ! Interior value 
    151          DO jj = 1, jpj 
    152             DO ji = 1, jpi 
    153                DO jk = mikt(ji,jj)+1, jpkm1 
    154                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    155                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
    156                   zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 
    157                END DO 
    158             END DO 
    159          END DO 
    160169 
    161170         ! total advective trend 
    … …  
    202211       
    203212         ! antidiffusive flux on k 
    204          zwz(:,:,1) = 0.e0         ! Surface value 
    205          ! 
    206          DO jj = 1, jpj 
    207             DO ji = 1, jpi 
    208                ik=mikt(ji,jj) 
    209                ! surface value 
    210                zwz(ji,jj,1:ik) = 0.e0 
    211                ! Interior value 
    212                DO jk = mikt(ji,jj)+1, jpkm1                     
     213         ! Interior value 
     214         DO jk = 2, jpkm1                     
     215            DO jj = 1, jpj 
     216               DO ji = 1, jpi 
    213217                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk) 
    214218               END DO 
    215219            END DO 
    216220         END DO 
     221         ! surface value 
     222         IF ( ln_isfcav ) THEN 
     223            DO jj = 1, jpj 
     224               DO ji = 1, jpi 
     225                  zwz(ji,jj,mikt(ji,jj)) = 0.e0 
     226               END DO 
     227            END DO 
     228         ELSE 
     229            zwz(:,:,1) = 0.e0 
     230         END IF 
    217231         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions 
    218232         CALL lbc_lnk( zwz, 'W',  1. ) 
    … …  
    358372 
    359373         ! upstream tracer flux in the k direction 
    360          ! Surface value 
    361          IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0._wp                        ! volume variable 
    362          ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface  
    363          ENDIF 
    364374         ! Interior value 
    365375         DO jk = 2, jpkm1 
    … …  
    372382            END DO 
    373383         END DO 
     384         ! Surface value 
     385         IF( lk_vvl ) THEN 
     386            IF ( ln_isfcav ) THEN 
     387               DO jj = 1, jpj 
     388                  DO ji = 1, jpi 
     389                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable +    isf 
     390                  END DO 
     391               END DO 
     392            ELSE 
     393               zwz(:,:,1) = 0.e0                              ! volume variable + no isf 
     394            END IF 
     395         ELSE 
     396            IF ( ln_isfcav ) THEN 
     397               DO jj = 1, jpj 
     398                  DO ji = 1, jpi 
     399                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface +    isf 
     400                  END DO 
     401               END DO 
     402            ELSE 
     403               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)                                               ! linear free surface + no isf 
     404            END IF 
     405         ENDIF 
    374406 
    375407         ! total advective trend 
    … …  
    580612         &        paft * tmask + zbig * ( 1._wp - tmask )  ) 
    581613 
    582       DO jj = 2, jpjm1 
    583          DO ji = fs_2, fs_jpim1   ! vector opt. 
    584             DO jk = mikt(ji,jj), jpkm1 
    585                ikm1 = MAX(jk-1,mikt(ji,jj)) 
    586                z2dtt = p2dt(jk) 
    587                 
     614      DO jk = 1, jpkm1 
     615         ikm1 = MAX(jk-1,1) 
     616         z2dtt = p2dt(jk) 
     617         DO jj = 2, jpjm1 
     618            DO ji = fs_2, fs_jpim1   ! vector opt. 
     619 
    588620               ! search maximum in neighbourhood 
    589621               zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r4990 r5120  
    290290      IF(lwp) WRITE(numout,*) '              homogeneous ocean T = ', zt0, ' S = ',zs0 
    291291 
     292      ! Initialisation of gtui/gtvi in case of no cavity 
     293      IF ( .NOT. ln_isfcav ) THEN 
     294         gtui(:,:,:) = 0.0_wp 
     295         gtvi(:,:,:) = 0.0_wp 
     296      END IF 
    292297      !                                        ! T & S profile (to be coded +namelist parameter 
    293298 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90

    r4990 r5120  
    116116               END DO 
    117117            END DO 
    118  
    119118            !                          !==  Laplacian  ==! 
    120119            ! 
    … …  
    125124               END DO 
    126125            END DO 
     126            ! 
    127127            IF( ln_zps ) THEN                ! set gradient at partial step level (last ocean level) 
    128128               DO jj = 1, jpjm1 
    … …  
    130130                     IF( mbku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn) 
    131131                     IF( mbkv(ji,jj) == jk )  ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn) 
    132                      ! (ISH) 
    133                      IF( miku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn) 
    134                      IF( mikv(ji,jj) == jk )  ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn) 
    135132                  END DO 
    136133               END DO 
    137134            ENDIF 
     135            ! (ISH) 
     136            IF( ln_zps .AND. ln_isfcav ) THEN ! set gradient at partial step level (first ocean level in a cavity) 
     137               DO jj = 1, jpjm1 
     138                  DO ji = 1, jpim1 
     139                     IF( miku(ji,jj) == MAX(jk,2) )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn) 
     140                     IF( mikv(ji,jj) == MAX(jk,2) )  ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn) 
     141                  END DO 
     142               END DO 
     143            ENDIF 
     144            ! 
    138145            DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient 
    139146               DO ji = fs_2, fs_jpim1   ! vector opt. 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r4990 r5120  
    106106      ! 
    107107      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     108      INTEGER  ::  ikt 
    108109      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars 
    109110      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
    … …  
    149150            END DO 
    150151         END DO 
     152 
     153         ! partial cell correction 
    151154         IF( ln_zps ) THEN      ! partial steps correction at the last ocean level  
    152155            DO jj = 1, jpjm1 
    153156               DO ji = 1, fs_jpim1   ! vector opt. 
    154157! IF useless if zpshde defines pgu everywhere 
    155                   IF (mbku(ji,jj) > 1) zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
    156                   IF (mbkv(ji,jj) > 1) zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    157                   ! (ISF) 
     158                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
     159                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
     160               END DO 
     161            END DO 
     162         ENDIF 
     163         IF( ln_zps .AND. ln_isfcav ) THEN      ! partial steps correction at the first wet level beneath a cavity 
     164            DO jj = 1, jpjm1 
     165               DO ji = 1, fs_jpim1   ! vector opt. 
    158166                  IF (miku(ji,jj) > 1) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)           
    159167                  IF (mikv(ji,jj) > 1) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)      
    160168               END DO 
    161169            END DO 
    162          ENDIF 
     170         END IF 
    163171 
    164172         !!---------------------------------------------------------------------- 
    165173         !!   II - horizontal trend  (full) 
    166174         !!---------------------------------------------------------------------- 
    167 !CDIR PARALLEL DO PRIVATE( zdk1t )  
    168          !                                                ! =============== 
    169          DO jj = 1, jpj                                 ! Horizontal slab 
    170             !                                             ! =============== 
    171             DO ji = 1, jpi   ! vector opt. 
    172                DO jk = mikt(ji,jj), jpkm1 
    173                ! 1. Vertical tracer gradient at level jk and jk+1 
    174                ! ------------------------------------------------ 
    175                ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
    176                   zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 
    177                ! 
    178                   IF( jk == mikt(ji,jj) ) THEN  ;   zdkt(ji,jj,jk) = zdk1t(ji,jj,jk) 
    179                   ELSE                          ;   zdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    180                   ENDIF 
     175!!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t )  
     176            ! 1. Vertical tracer gradient at level jk and jk+1 
     177            ! ------------------------------------------------ 
     178         !  
     179         ! interior value  
     180         DO jk = 2, jpkm1                
     181            DO jj = 1, jpj 
     182               DO ji = 1, jpi   ! vector opt. 
     183                  zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn  ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) 
     184                  ! 
     185                  zdkt(ji,jj,jk)  = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn  ) ) * wmask(ji,jj,jk) 
    181186               END DO 
    182187            END DO 
    183188         END DO 
    184  
    185             ! 2. Horizontal fluxes 
    186             ! --------------------    
    187          DO jj = 1 , jpjm1 
    188             DO ji = 1, fs_jpim1   ! vector opt. 
    189                DO jk = mikt(ji,jj), jpkm1 
     189         ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
     190         zdk1t(:,:,1) = ( ptb(:,:,1,jn  ) - ptb(:,:,2,jn) ) * wmask(:,:,2) 
     191         zdkt (:,:,1) = zdk1t(:,:,1) 
     192         IF ( ln_isfcav ) THEN 
     193            DO jj = 1, jpj 
     194               DO ji = 1, jpi   ! vector opt. 
     195                  ikt = mikt(ji,jj) ! surface level 
     196                  zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 
     197                  zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 
     198               END DO 
     199            END DO 
     200         END IF 
     201 
     202         ! 2. Horizontal fluxes 
     203         ! --------------------    
     204         DO jk = 1, jpkm1 
     205            DO jj = 1 , jpjm1 
     206               DO ji = 1, fs_jpim1   ! vector opt. 
    190207                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
    191208                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
    … …  
    208225               END DO 
    209226            END DO 
    210          END DO 
    211227 
    212228            ! II.4 Second derivative (divergence) and add to the general trend 
    213229            ! ---------------------------------------------------------------- 
    214          DO jj = 2 , jpjm1 
    215             DO ji = fs_2, fs_jpim1   ! vector opt. 
    216                DO jk = mikt(ji,jj), jpkm1 
    217                   zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     230            DO jj = 2 , jpjm1 
     231               DO ji = fs_2, fs_jpim1   ! vector opt. 
     232                  zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    218233                  ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  ) 
    219234                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
    … …  
    278293            DO jj = 2, jpjm1 
    279294               DO ji = fs_2, fs_jpim1   ! vector opt. 
    280                   zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     295                  zcoef0 = - fsahtw(ji,jj,jk) * wmask(ji,jj,jk) 
    281296                  ! 
    282297                  zmsku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      & 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90

    r4990 r5120  
    102102               END DO 
    103103            END DO 
    104             IF( ln_zps ) THEN      ! set gradient at partial step level 
     104            IF( ln_zps ) THEN      ! set gradient at partial step level for the last ocean cell 
    105105               DO jj = 1, jpjm1 
    106106                  DO ji = 1, fs_jpim1   ! vector opt. 
    … …  
    116116                        ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn) 
    117117                     ENDIF 
    118                       
    119                      ! (ISH) 
     118                  END DO 
     119               END DO 
     120            ENDIF 
     121            ! (ISH) 
     122            IF( ln_zps .AND. ln_isfcav ) THEN      ! set gradient at partial step level for the first ocean cell 
     123                                                   ! into a cavity 
     124               DO jj = 1, jpjm1 
     125                  DO ji = 1, fs_jpim1   ! vector opt. 
    120126                     ! ice shelf level level MAX(2,jk) => only where ice shelf 
    121127                     iku = miku(ji,jj)  

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r4990 r5120  
    99   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
    1010   !!             -   !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC 
     11   !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing  
    1112   !!---------------------------------------------------------------------- 
    1213 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90

    r4990 r5120  
    122122            DO jj=1, jpj 
    123123               DO ji=1, jpi 
    124                   zwt(ji,jj,1:mikt(ji,jj)) = 0._wp 
     124                  zwt(ji,jj,1) = 0._wp 
    125125               END DO 
    126126            END DO 
    … …  
    184184            DO jj = 2, jpjm1 
    185185               DO ji = fs_2, fs_jpim1 
    186                   zwt(ji,jj,1:mikt(ji,jj)) = zwd(ji,jj,1:mikt(ji,jj)) 
    187                   DO jk = mikt(ji,jj)+1, jpkm1 
     186                  zwt(ji,jj,1) = zwd(ji,jj,1) 
     187               END DO 
     188            END DO 
     189            DO jk = 2, jpkm1 
     190               DO jj = 2, jpjm1 
     191                  DO ji = fs_2, fs_jpim1 
    188192                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
    189193                  END DO 
    … …  
    196200         DO jj = 2, jpjm1 
    197201            DO ji = fs_2, fs_jpim1 
    198                ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,mikt(ji,jj)) 
    199                ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,mikt(ji,jj)) 
    200                pta(ji,jj,mikt(ji,jj),jn) = ze3tb * ptb(ji,jj,mikt(ji,jj),jn)                     & 
    201                   &                      + p2dt(mikt(ji,jj)) * ze3tn * pta(ji,jj,mikt(ji,jj),jn) 
    202                DO jk = mikt(ji,jj)+1, jpkm1 
     202               ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 
     203               ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1) 
     204               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn)                     & 
     205                  &                      + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 
     206            END DO 
     207         END DO 
     208         DO jk = 2, jpkm1 
     209            DO jj = 2, jpjm1 
     210               DO ji = fs_2, fs_jpim1 
    203211                  ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk) 
    204212                  ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t  (ji,jj,jk) 
    … …  
    213221            DO ji = fs_2, fs_jpim1 
    214222               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
    215                DO jk = jpk-2, mikt(ji,jj), -1 
     223            END DO 
     224         END DO 
     225         DO jk = jpk-2, 1, -1 
     226            DO jj = 2, jpjm1 
     227               DO ji = fs_2, fs_jpim1 
    216228                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   & 
    217229                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r4990 r5120  
    88   !!             -   !  2004-03  (C. Ethe)  adapted for passive tracers 
    99   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA  
     10   !!            3.6  !  2014-11  (P. Mathiot) Add zps_hde_isf (needed to open a cavity) 
    1011   !!====================================================================== 
    1112    
    … …  
    2728   PRIVATE 
    2829 
    29    PUBLIC   zps_hde    ! routine called by step.F90 
     30   PUBLIC   zps_hde     ! routine called by step.F90 
     31   PUBLIC   zps_hde_isf ! routine called by step.F90 
    3032 
    3133   !! * Substitutions 
    … …  
    4042 
    4143   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   & 
     44      &                          prd, pgru, pgrv    ) 
     45      !!---------------------------------------------------------------------- 
     46      !!                     ***  ROUTINE zps_hde  *** 
     47      !!                     
     48      !! ** Purpose :   Compute the horizontal derivative of T, S and rho 
     49      !!      at u- and v-points with a linear interpolation for z-coordinate 
     50      !!      with partial steps. 
     51      !! 
     52      !! ** Method  :   In z-coord with partial steps, scale factors on last  
     53      !!      levels are different for each grid point, so that T, S and rd  
     54      !!      points are not at the same depth as in z-coord. To have horizontal 
     55      !!      gradients again, we interpolate T and S at the good depth :  
     56      !!      Linear interpolation of T, S    
     57      !!         Computation of di(tb) and dj(tb) by vertical interpolation: 
     58      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~ 
     59      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~ 
     60      !!         This formulation computes the two cases: 
     61      !!                 CASE 1                   CASE 2   
     62      !!         k-1  ___ ___________   k-1   ___ ___________ 
     63      !!                    Ti  T~                  T~  Ti+1 
     64      !!                  _____                        _____ 
     65      !!         k        |   |Ti+1     k           Ti |   | 
     66      !!                  |   |____                ____|   | 
     67      !!              ___ |   |   |           ___  |   |   | 
     68      !!                   
     69      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then 
     70      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1) 
     71      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  ) 
     72      !!          or 
     73      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then 
     74      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i ) 
     75      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) ) 
     76      !!          Idem for di(s) and dj(s)           
     77      !! 
     78      !!      For rho, we call eos which will compute rd~(t~,s~) at the right 
     79      !!      depth zh from interpolated T and S for the different formulations 
     80      !!      of the equation of state (eos). 
     81      !!      Gradient formulation for rho : 
     82      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~ 
     83      !! 
     84      !! ** Action  : compute for top interfaces 
     85      !!              - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 
     86      !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 
     87      !!---------------------------------------------------------------------- 
     88      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
     89      INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers 
     90      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields 
     91      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
     92      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
     93      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom) 
     94      ! 
     95      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
     96      INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points 
     97      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars 
     98      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
     99      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !  
     100      !!---------------------------------------------------------------------- 
     101      ! 
     102      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde') 
     103      ! 
     104      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
     105      zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
     106      zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
     107      ! 
     108      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
     109         ! 
     110         DO jj = 1, jpjm1 
     111            DO ji = 1, jpim1 
     112               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
     113               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
     114               ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
     115               ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     116               ! 
     117               ! i- direction 
     118               IF( ze3wu >= 0._wp ) THEN      ! case 1 
     119                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
     120                  ! interpolated values of tracers 
     121                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     122                  ! gradient of  tracers 
     123                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     124               ELSE                           ! case 2 
     125                  zmaxu = -ze3wu / fse3w(ji,jj,iku) 
     126                  ! interpolated values of tracers 
     127                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     128                  ! gradient of tracers 
     129                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     130               ENDIF 
     131               ! 
     132               ! j- direction 
     133               IF( ze3wv >= 0._wp ) THEN      ! case 1 
     134                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
     135                  ! interpolated values of tracers 
     136                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     137                  ! gradient of tracers 
     138                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     139               ELSE                           ! case 2 
     140                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
     141                  ! interpolated values of tracers 
     142                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     143                  ! gradient of tracers 
     144                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     145               ENDIF 
     146            END DO 
     147         END DO 
     148         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     149         ! 
     150      END DO 
     151 
     152      ! horizontal derivative of density anomalies (rd) 
     153      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
     154         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
     155         DO jj = 1, jpjm1 
     156            DO ji = 1, jpim1 
     157               iku = mbku(ji,jj) 
     158               ikv = mbkv(ji,jj) 
     159               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
     160               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     161               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1 
     162               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2 
     163               ENDIF 
     164               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1 
     165               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2 
     166               ENDIF 
     167            END DO 
     168         END DO 
     169 
     170         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
     171         ! step and store it in  zri, zrj for each  case 
     172         CALL eos( zti, zhi, zri )   
     173         CALL eos( ztj, zhj, zrj ) 
     174 
     175         ! Gradient of density at the last level  
     176         DO jj = 1, jpjm1 
     177            DO ji = 1, jpim1 
     178               iku = mbku(ji,jj) 
     179               ikv = mbkv(ji,jj) 
     180               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
     181               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     182               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
     183               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2 
     184               ENDIF 
     185               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1 
     186               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2 
     187               ENDIF 
     188            END DO 
     189         END DO 
     190         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
     191         ! 
     192      END IF 
     193      ! 
     194      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde') 
     195      ! 
     196   END SUBROUTINE zps_hde 
     197   ! 
     198   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv,   & 
    42199      &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  & 
    43       &                   sgtu, sgtv, sgru, sgrv, smru, smrv, sgzu, sgzv, sge3ru, sge3rv ) 
     200      &                   pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 
    44201      !!---------------------------------------------------------------------- 
    45202      !!                     ***  ROUTINE zps_hde  *** 
    … …  
    82239      !! 
    83240      !! ** Action  : compute for top and bottom interfaces 
    84       !!              - pgtu, pgtv, sgtu, sgtv: horizontal gradient of tracer at u- & v-points 
    85       !!              - pgru, pgrv, sgru, sgtv: horizontal gradient of rho (if present) at u- & v-points 
    86       !!              - pmru, pmrv, smru, smrv: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 
    87       !!              - pgzu, pgzv, sgzu, sgzv: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 
    88       !!              - pge3ru, pge3rv, sge3ru, sge3rv: horizontal gradient of rho weighted by local e3w at u- & v-points  
     241      !!              - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points 
     242      !!              - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 
     243      !!              - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 
     244      !!              - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 
     245      !!              - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points  
    89246      !!---------------------------------------------------------------------- 
    90247      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
    … …  
    92249      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields 
    93250      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
    94       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  sgtu, sgtv  ! hor. grad. of stra at u- & v-pts (ISF) 
     251      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi  ! hor. grad. of stra at u- & v-pts (ISF) 
    95252      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    96253      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom) 
    … …  
    98255      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom) 
    99256      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom) 
    100       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sgru, sgrv      ! hor. grad of prd at u- & v-pts (top) 
    101       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  smru, smrv      ! hor. sum  of prd at u- & v-pts (top) 
    102       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sgzu, sgzv      ! hor. grad of z   at u- & v-pts (top) 
    103       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sge3ru, sge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 
     257      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi      ! hor. grad of prd at u- & v-pts (top) 
     258      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui, pmrvi      ! hor. sum  of prd at u- & v-pts (top) 
     259      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui, pgzvi      ! hor. grad of z   at u- & v-pts (top) 
     260      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui, pge3rvi  ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 
    104261      ! 
    105262      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    … …  
    110267      !!---------------------------------------------------------------------- 
    111268      ! 
    112       IF( nn_timing == 1 )  CALL timing_start( 'zps_hde') 
     269      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf') 
    113270      ! 
    114271      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
    115       sgtu(:,:,:)=0.0_wp ; sgtv(:,:,:)=0.0_wp ; 
     272      pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ; 
    116273      zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
    117274      zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
    … …  
    256413                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 
    257414                  ! gradient of tracers 
    258                   sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     415                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    259416               ELSE                           ! case 2 
    260417                  zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 
    … …  
    262419                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 
    263420                  ! gradient of  tracers 
    264                   sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     421                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    265422               ENDIF 
    266423               ! 
    … …  
    271428                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 
    272429                  ! gradient of tracers 
    273                   sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     430                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    274431               ELSE                           ! case 2 
    275432                  zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1) 
    … …  
    277434                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 
    278435                  ! gradient of tracers 
    279                   sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     436                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    280437               ENDIF 
    281438            END DO!! 
    282439         END DO!! 
    283          CALL lbc_lnk( sgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( sgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     440         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    284441         ! 
    285442      END DO 
    … …  
    287444      ! horizontal derivative of density anomalies (rd) 
    288445      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    289          sgru(:,:)  =0.0_wp ; sgrv(:,:)  =0.0_wp ; 
    290          sgzu(:,:)  =0.0_wp ; sgzv(:,:)  =0.0_wp ; 
    291          smru(:,:)  =0.0_wp ; smru(:,:)  =0.0_wp ; 
    292          sge3ru(:,:)=0.0_wp ; sge3rv(:,:)=0.0_wp ; 
     446         pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ; 
     447         pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ; 
     448         pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ; 
     449         pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 
    293450 
    294451         DO jj = 1, jpjm1 
    … …  
    321478               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 
    322479               IF( ze3wu >= 0._wp ) THEN 
    323                  sgzu  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 
    324                  sgru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1 
    325                  smru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1  
    326                  sge3ru(ji,jj) = umask(ji,jj,iku+1)                                                                  & 
     480                 pgzui  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 
     481                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1 
     482                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1  
     483                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  & 
    327484                                * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   & 
    328485                                   - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1 
    329486               ELSE 
    330                  sgzu  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 
    331                  sgru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2 
    332                  smru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2 
    333                  sge3ru(ji,jj) = umask(ji,jj,iku+1)                                                                   & 
     487                 pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 
     488                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2 
     489                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2 
     490                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   & 
    334491                                * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  & 
    335492                                   -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2 
    336493               ENDIF 
    337494               IF( ze3wv >= 0._wp ) THEN 
    338                  sgzv  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)  
    339                  sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1 
    340                  smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1 
    341                  sge3rv(ji,jj) = vmask(ji,jj,ikv+1)                                                                  &  
     495                 pgzvi  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)  
     496                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1 
     497                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1 
     498                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  &  
    342499                                * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  & 
    343500                                   - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1 
    344501                                  ! + 2 due to the formulation in density and not in anomalie in hpg sco 
    345502               ELSE 
    346                  sgzv  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 
    347                  sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2 
    348                  smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2 
    349                  sge3rv(ji,jj) = vmask(ji,jj,ikv+1)                                                                   & 
     503                 pgzvi  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 
     504                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2 
     505                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2 
     506                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   & 
    350507                                * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 
    351508                                   -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2 
    … …  
    353510            END DO 
    354511         END DO 
    355          CALL lbc_lnk( sgru   , 'U', -1. )   ;   CALL lbc_lnk( sgrv   , 'V', -1. )   ! Lateral boundary conditions 
    356          CALL lbc_lnk( smru   , 'U',  1. )   ;   CALL lbc_lnk( smrv   , 'V',  1. )   ! Lateral boundary conditions 
    357          CALL lbc_lnk( sgzu   , 'U', -1. )   ;   CALL lbc_lnk( sgzv   , 'V', -1. )   ! Lateral boundary conditions 
    358          CALL lbc_lnk( sge3ru , 'U', -1. )   ;   CALL lbc_lnk( sge3rv , 'V', -1. )   ! Lateral boundary conditions 
     512         CALL lbc_lnk( pgrui   , 'U', -1. )   ;   CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions 
     513         CALL lbc_lnk( pmrui   , 'U',  1. )   ;   CALL lbc_lnk( pmrvi   , 'V',  1. )   ! Lateral boundary conditions 
     514         CALL lbc_lnk( pgzui   , 'U', -1. )   ;   CALL lbc_lnk( pgzvi   , 'V', -1. )   ! Lateral boundary conditions 
     515         CALL lbc_lnk( pge3rui , 'U', -1. )   ;   CALL lbc_lnk( pge3rvi , 'V', -1. )   ! Lateral boundary conditions 
    359516         ! 
    360517      END IF   
    361518      ! 
    362       IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde') 
    363       ! 
    364    END SUBROUTINE zps_hde 
    365  
     519      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde_isf') 
     520      ! 
     521   END SUBROUTINE zps_hde_isf 
    366522   !!====================================================================== 
    367523END MODULE zpshde 

  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r4990 r5120  
    120120                  zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 
    121121                  zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max) 
    122 ! (ISF) 
    123                   ikbt = mikt(ji,jj) 
    124 ! JC: possible WAD implementation should modify line below if layers vanish 
    125                   ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
    126                   ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 
    127                   ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max) 
    128  
    129122               END DO 
    130123            END DO 
     124! (ISF) 
     125            IF ( ln_isfcav ) THEN 
     126               DO jj = 1, jpj 
     127                  DO ji = 1, jpi 
     128                     ikbt = mikt(ji,jj) 
     129! JC: possible WAD implementation should modify line below if layers vanish 
     130                     ztmp = (1-tmask(ji,jj,1)) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
     131                     ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 
     132                     ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max) 
     133                  END DO 
     134               END DO 
     135            END IF 
    131136         !    
    132137         ELSE 
    … …  
    152157               ! 
    153158               ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
    154                IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 
    155                   bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) )   & 
    156                                &            + ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) ) & 
    157                                &          * zecu * (1._wp - umask(ji,jj,1)) 
    158                END IF 
    159                IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 
    160                   bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) )   & 
    161                                &            + ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) ) & 
    162                                &          * zecv * (1._wp - vmask(ji,jj,1)) 
    163                END IF 
    164                ! (ISF) ======================================================================== 
    165                ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points  
    166                ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points) 
    167                ! 
    168                zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     & 
    169                   &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  ) 
    170                zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     & 
    171                   &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  ) 
    172                ! 
    173                zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 ) 
    174                zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 
    175                ! 
    176                tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1)) 
    177                tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1)) 
    178                ! (ISF) END ==================================================================== 
    179                ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
    180                IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 
    181                   tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) )   & 
    182                                &            + ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) ) & 
    183                                &          * zecu * (1._wp - umask(ji,jj,1)) 
    184                END IF 
    185                IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 
    186                   tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) )   & 
    187                                &            + ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) ) & 
    188                                &          * zecv * (1._wp - vmask(ji,jj,1)) 
     159               IF ( ln_isfcav ) THEN 
     160                  IF ( miku(ji,jj) + 1 .GE. mbku(ji,jj) ) THEN 
     161                     bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) )   & 
     162                                  &            + ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) ) & 
     163                                  &          * zecu * (1._wp - umask(ji,jj,1)) 
     164                  END IF 
     165                  IF ( mikv(ji,jj) + 1 .GE. mbkv(ji,jj) ) THEN 
     166                     bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) )   & 
     167                                  &            + ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) ) & 
     168                                  &          * zecv * (1._wp - vmask(ji,jj,1)) 
     169                  END IF 
    189170               END IF 
    190171            END DO 
    191172         END DO 
     173         IF ( ln_isfcav ) THEN 
     174            DO jj = 2, jpjm1 
     175               DO ji = 2, jpim1 
     176                  ! (ISF) ======================================================================== 
     177                  ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points  
     178                  ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points) 
     179                  ! 
     180                  zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     & 
     181                     &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  ) 
     182                  zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     & 
     183                     &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  ) 
     184              ! 
     185                  zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 ) 
     186                  zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 
     187              ! 
     188                  tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1)) 
     189                  tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1)) 
     190              ! (ISF) END ==================================================================== 
     191              ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
     192                  IF ( miku(ji,jj) + 1 .GE. mbku(ji,jj) ) THEN 
     193                     tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) )   & 
     194                                  &            + ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) ) & 
     195                                  &          * zecu * (1._wp - umask(ji,jj,1)) 
     196                  END IF 
     197                  IF ( mikv(ji,jj) + 1 .GE. mbkv(ji,jj) ) THEN 
     198                     tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) )   & 
     199                                  &            + ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) ) & 
     200                                  &          * zecv * (1._wp - vmask(ji,jj,1)) 
     201                  END IF 
     202               END DO 
     203            END DO 
     204         END IF 
    192205         ! 
    193206         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition 

  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r4990 r5120  
    156156         END DO 
    157157         ! mask zmsk in order to have avt and avs masked 
    158          zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk) 
     158         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk) 
    159159 
    160160 
    … …  
    191191               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    & 
    192192                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   & 
    193                   &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk) 
     193                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * wumask(ji,jj,jk) 
    194194               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    & 
    195195                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   & 
    196                   &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * vmask(ji,jj,jk) 
     196                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * wvmask(ji,jj,jk) 
    197197            END DO 
    198198         END DO 
    … …  
    255255      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' ) 
    256256      !                               ! initialization to masked Kz 
    257       avs(:,:,:) = rn_avt0 * tmask(:,:,:)  
     257      avs(:,:,:) = rn_avt0 * wmask(:,:,:)  
    258258      ! 
    259259   END SUBROUTINE zdf_ddm_init 

  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

    r4990 r5120  
    1414   !!---------------------------------------------------------------------- 
    1515   USE par_oce         ! mesh and scale factors 
    16    USE sbc_oce         ! surface module (only for nn_isf in the option compatibility test) 
    1716   USE ldftra_oce      ! ocean active tracers: lateral physics 
    1817   USE ldfdyn_oce      ! ocean dynamics lateral physics 
    … …  
    118117      IF( ioptio == 0 .OR. ioptio > 1 .AND. .NOT. lk_esopa )   & 
    119118         &   CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' ) 
    120       IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. nn_isf .NE. 0 )   & 
     119      IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. ln_isfcav )   & 
    121120         &   CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' ) 
    122121      ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r5112 r5120  
    2626   !!                 !                                + cleaning of the parameters + bugs correction 
    2727   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     28   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    2829   !!---------------------------------------------------------------------- 
    2930#if defined key_zdftke   ||   defined key_esopa 
    … …  
    236237      zfact3 = 0.5_wp       * rn_ediss 
    237238      ! 
     239      ! 
    238240      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    239241      !                     !  Surface boundary condition on tke 
    240242      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     243      IF ( ln_isfcav ) THEN 
     244         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin 
     245            DO ji = fs_2, fs_jpim1   ! vector opt. 
     246               en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1) 
     247            END DO 
     248         END DO 
     249      END IF 
    241250      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    242251         DO ji = fs_2, fs_jpim1   ! vector opt. 
    243             IF (mikt(ji,jj) .GT. 1) THEN 
    244                en(ji,jj,mikt(ji,jj))=rn_emin 
    245             ELSE 
    246                en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    247             END IF 
     252            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    248253         END DO 
    249254      END DO 
    … …  
    301306         END DO 
    302307         zcof = 0.016 / SQRT( zrhoa * zcdrag ) 
     308!CDIR NOVERRCHK 
    303309         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en 
    304             DO jj = 2, jpjm1 
     310!CDIR NOVERRCHK 
     311            DO jj = 2, jpjm1 
     312!CDIR NOVERRCHK 
    305313               DO ji = fs_2, fs_jpim1   ! vector opt. 
    306314                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
    … …  
    309317                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    310318                  !                                           ! TKE Langmuir circulation source term 
    311                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk) 
     319                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    312320               END DO 
    313321            END DO 
    … …  
    328336               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   & 
    329337                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &  
    330                   &           / (  fse3uw_n(ji,jj,jk)         & 
    331                   &              * fse3uw_b(ji,jj,jk) ) 
     338                  &                            / (  fse3uw_n(ji,jj,jk)               & 
     339                  &                              *  fse3uw_b(ji,jj,jk) ) 
    332340               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   & 
    333341                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   & 
    … …  
    338346      END DO 
    339347      ! 
    340       DO jj = 2, jpjm1 
    341          DO ji = fs_2, fs_jpim1   ! vector opt. 
    342             DO jk = mikt(ji,jj)+1, jpkm1           !* Matrix and right hand side in en 
     348      DO jk = 2, jpkm1           !* Matrix and right hand side in en 
     349         DO jj = 2, jpjm1 
     350            DO ji = fs_2, fs_jpim1   ! vector opt. 
    343351               zcof   = zfact1 * tmask(ji,jj,jk) 
    344352               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal 
    … …  
    357365               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
    358366                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
    359                   &                                 * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
    360             END DO 
    361             !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
    362             DO jk = mikt(ji,jj)+2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
     367                  &                                 * wmask(ji,jj,jk) 
     368            END DO 
     369         END DO 
     370      END DO 
     371      !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
     372      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
     373         DO jj = 2, jpjm1 
     374            DO ji = fs_2, fs_jpim1    ! vector opt. 
    363375               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    364376            END DO 
    365             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    366             zd_lw(ji,jj,mikt(ji,jj)+1) = en(ji,jj,mikt(ji,jj)+1) - zd_lw(ji,jj,mikt(ji,jj)+1) * en(ji,jj,mikt(ji,jj))    ! Surface boudary conditions on tke 
    367             ! 
    368             DO jk = mikt(ji,jj)+2, jpkm1 
     377         END DO 
     378      END DO 
     379      ! 
     380      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
     381      DO jj = 2, jpjm1 
     382         DO ji = fs_2, fs_jpim1   ! vector opt. 
     383            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
     384         END DO 
     385      END DO 
     386      DO jk = 3, jpkm1 
     387         DO jj = 2, jpjm1 
     388            DO ji = fs_2, fs_jpim1    ! vector opt. 
    369389               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 
    370390            END DO 
    371             ! 
    372             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
     391         END DO 
     392      END DO 
     393      ! 
     394      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
     395      DO jj = 2, jpjm1 
     396         DO ji = fs_2, fs_jpim1   ! vector opt. 
    373397            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
    374             ! 
    375             DO jk = jpk-2, mikt(ji,jj)+1, -1 
     398         END DO 
     399      END DO 
     400      DO jk = jpk-2, 2, -1 
     401         DO jj = 2, jpjm1 
     402            DO ji = fs_2, fs_jpim1    ! vector opt. 
    376403               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    377404            END DO 
    378             ! 
    379             DO jk = mikt(ji,jj), jpkm1                             ! set the minimum value of tke 
    380                en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 
     405         END DO 
     406      END DO 
     407      DO jk = 2, jpkm1                             ! set the minimum value of tke 
     408         DO jj = 2, jpjm1 
     409            DO ji = fs_2, fs_jpim1   ! vector opt. 
     410               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 
    381411            END DO 
    382412         END DO 
    … …  
    391421               DO ji = fs_2, fs_jpim1   ! vector opt. 
    392422                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    393                      &                                 * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
     423                     &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    394424               END DO 
    395425            END DO 
    … …  
    400430               jk = nmln(ji,jj) 
    401431               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    402                   &                                 * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
     432                  &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    403433            END DO 
    404434         END DO 
    … …  
    416446                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    417447                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    418                      &                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1) 
     448                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    419449               END DO 
    420450            END DO 
    … …  
    484514      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2) 
    485515      ! 
     516      ! initialisation of interior minimum value (avoid a 2d loop with mikt) 
     517      zmxlm(:,:,:)  = rmxl_min     
     518      zmxld(:,:,:)  = rmxl_min 
     519      ! 
    486520      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    487521         DO jj = 2, jpjm1 
    488522            DO ji = fs_2, fs_jpim1 
    489                IF (mikt(ji,jj) .GT. 1) THEN 
    490                   zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min 
    491                ELSE 
    492                   zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
    493                   zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) ) 
    494                END IF 
     523               zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
     524               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
    495525            END DO 
    496526         END DO 
    497527      ELSE  
    498          DO jj = 2, jpjm1 
    499             DO ji = fs_2, fs_jpim1                         ! surface set to the minimum value 
    500                zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) 
    501             END DO 
    502          END DO 
     528         zmxlm(:,:,1) = rn_mxl0 
    503529      ENDIF 
    504       zmxlm(:,:,jpk)  = rmxl_min     ! last level set to the interior minium value 
    505       ! 
    506 !CDIR NOVERRCHK 
    507       DO jj = 2, jpjm1 
    508 !CDIR NOVERRCHK 
    509          DO ji = fs_2, fs_jpim1   ! vector opt. 
    510             !CDIR NOVERRCHK 
    511             DO jk = mikt(ji,jj)+1, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
     530      ! 
     531!CDIR NOVERRCHK 
     532      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
     533!CDIR NOVERRCHK 
     534         DO jj = 2, jpjm1 
     535!CDIR NOVERRCHK 
     536            DO ji = fs_2, fs_jpim1   ! vector opt. 
    512537               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    513                zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  ) 
    514             END DO 
    515             zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj))   ! surface set to the minimum value  
     538               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 
     539            END DO 
    516540         END DO 
    517541      END DO 
    … …  
    519543      !                     !* Physical limits for the mixing length 
    520544      ! 
    521       zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the zmxlm   value 
     545      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value  
    522546      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value 
    523547      ! 
    524548      SELECT CASE ( nn_mxl ) 
    525549      ! 
     550      ! where wmask = 0 set zmxlm == fse3w 
    526551      CASE ( 0 )           ! bounded by the distance to surface and bottom 
    527          DO jj = 2, jpjm1 
    528             DO ji = fs_2, fs_jpim1   ! vector opt. 
    529                DO jk = mikt(ji,jj)+1, jpkm1 
     552         DO jk = 2, jpkm1 
     553            DO jj = 2, jpjm1 
     554               DO ji = fs_2, fs_jpim1   ! vector opt. 
    530555                  zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   & 
    531556                  &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 
    532                   zmxlm(ji,jj,jk) = zemxl 
    533                   zmxld(ji,jj,jk) = zemxl 
     557                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 
     558                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 
     559                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 
    534560               END DO 
    535561            END DO 
    … …  
    537563         ! 
    538564      CASE ( 1 )           ! bounded by the vertical scale factor 
    539          DO jj = 2, jpjm1 
    540             DO ji = fs_2, fs_jpim1   ! vector opt. 
    541                DO jk = mikt(ji,jj)+1, jpkm1 
     565         DO jk = 2, jpkm1 
     566            DO jj = 2, jpjm1 
     567               DO ji = fs_2, fs_jpim1   ! vector opt. 
    542568                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 
    543569                  zmxlm(ji,jj,jk) = zemxl 
    … …  
    548574         ! 
    549575      CASE ( 2 )           ! |dk[xml]| bounded by e3t : 
    550          DO jj = 2, jpjm1 
    551             DO ji = fs_2, fs_jpim1   ! vector opt. 
    552                DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : 
     576         DO jk = 2, jpkm1         ! from the surface to the bottom : 
     577            DO jj = 2, jpjm1 
     578               DO ji = fs_2, fs_jpim1   ! vector opt. 
    553579                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    554580               END DO 
    555                DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : 
     581            END DO 
     582         END DO 
     583         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : 
     584            DO jj = 2, jpjm1 
     585               DO ji = fs_2, fs_jpim1   ! vector opt. 
    556586                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    557587                  zmxlm(ji,jj,jk) = zemxl 
    … …  
    562592         ! 
    563593      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t : 
    564          DO jj = 2, jpjm1 
    565             DO ji = fs_2, fs_jpim1   ! vector opt. 
    566                DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : lup 
     594         DO jk = 2, jpkm1         ! from the surface to the bottom : lup 
     595            DO jj = 2, jpjm1 
     596               DO ji = fs_2, fs_jpim1   ! vector opt. 
    567597                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    568598               END DO 
    569                DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : ldown 
     599            END DO 
     600         END DO 
     601         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown 
     602            DO jj = 2, jpjm1 
     603               DO ji = fs_2, fs_jpim1   ! vector opt. 
    570604                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    571605               END DO 
    … …  
    604638               zsqen = SQRT( en(ji,jj,jk) ) 
    605639               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
    606                avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk) 
    607                avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
     640               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
     641               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
    608642               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 
    609643            END DO 
    … …  
    612646      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    613647      ! 
    614       DO jj = 2, jpjm1 
    615          DO ji = fs_2, fs_jpim1   ! vector opt. 
    616             DO jk = miku(ji,jj)+1, jpkm1            !* vertical eddy viscosity at u- and v-points 
    617                avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
    618             END DO 
    619             DO jk = mikv(ji,jj)+1, jpkm1 
    620                avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
     648      DO jk = 2, jpkm1            !* vertical eddy viscosity at wu- and wv-points 
     649         DO jj = 2, jpjm1 
     650            DO ji = fs_2, fs_jpim1   ! vector opt. 
     651               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
     652               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
    621653            END DO 
    622654         END DO 
    … …  
    625657      ! 
    626658      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
    627          DO jj = 2, jpjm1 
    628             DO ji = fs_2, fs_jpim1   ! vector opt. 
    629                DO jk = mikt(ji,jj)+1, jpkm1 
     659         DO jk = 2, jpkm1 
     660            DO jj = 2, jpjm1 
     661               DO ji = fs_2, fs_jpim1   ! vector opt. 
    630662                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 
    631663                  !                                          ! shear 
    … …  
    639671!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!) 
    640672!!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  ) 
    641                   avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
     673                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
    642674# if defined key_c1d 
    643                   e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number 
    644                   e_ric(ji,jj,jk) = zri   * tmask(ji,jj,jk)  ! c1d config. : save Ri 
     675                  e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number 
     676                  e_ric(ji,jj,jk) = zri   * wmask(ji,jj,jk)  ! c1d config. : save Ri 
    645677# endif 
    646678              END DO 
    … …  
    749781      !                               !* set vertical eddy coef. to the background value 
    750782      DO jk = 1, jpk 
    751          avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 
    752          avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 
    753          avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 
    754          avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 
     783         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
     784         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
     785         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 
     786         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 
    755787      END DO 
    756788      dissl(:,:,:) = 1.e-12_wp 
    … …  
    814846           en(:,:,:) = rn_emin * tmask(:,:,:) 
    815847           DO jk = 1, jpk                           ! set the Kz to the background value 
    816               avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 
    817               avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 
    818               avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 
    819               avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 
     848              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
     849              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
     850              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 
     851              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 
    820852           END DO 
    821853        ENDIF 

  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r5021 r5120  
    126126      zkz(:,:) = 0.e0               !* Associated potential energy consummed over the whole water column 
    127127      DO jk = 2, jpkm1 
    128          zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk) * tmask(:,:,jk-1) 
     128         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk) * wmask(:,:,jk) 
    129129      END DO 
    130130 
    … …  
    135135      END DO 
    136136 
    137       DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
    138          DO ji = 1, jpi 
    139             DO jk = mikt(ji,jj)+1, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 
    140                zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s 
     137      DO jk = 2, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 
     138         DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
     139            DO ji = 1, jpi 
     140               zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk)  !kz max = 300 cm2/s 
    141141            END DO 
    142142         END DO 
    … …  
    166166      !                          !   Update  mixing coefs  !                           
    167167      !                          ! ----------------------- ! 
    168       DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
    169          DO ji = 1, jpi 
    170             DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
    171                avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) 
    172                avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) 
     168      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
     169         DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
     170            DO ji = 1, jpi 
     171               avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk) 
     172               avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk) 
    173173            END DO 
    174174         END DO 
    175175      END DO 
    176176       
    177       DO jj = 2, jpjm1 
    178          DO ji = fs_2, fs_jpim1  ! vector opt. 
    179             DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
    180                avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
    181                avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
     177      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
     178         DO jj = 2, jpjm1 
     179            DO ji = fs_2, fs_jpim1  ! vector opt. 
     180               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
     181               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
    182182            END DO 
    183183         END DO 
    … …  
    457457         ztpc = 0.e0 
    458458         zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:) 
    459          DO jj = 1, jpj 
    460             DO ji = 1, jpi 
    461                DO jk= mikt(ji,jj)+1, jpkm1 
    462                   ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 
     459         DO jk= 2, jpkm1 
     460            DO jj = 1, jpj 
     461               DO ji = 1, jpi 
     462                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
    463463               END DO 
    464464            END DO 
    … …  
    473473         zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )    
    474474         zkz(:,:) = 0.e0 
    475          DO jj = 1, jpj 
    476             DO ji = 1, jpi 
    477                DO jk = mikt(ji,jj)+1, jpkm1 
    478                zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* tmask(ji,jj,jk) 
     475         DO jk = 2, jpkm1 
     476            DO jj = 1, jpj 
     477               DO ji = 1, jpi 
     478                  zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* wmask(ji,jj,jk) 
    479479               END DO 
    480480            END DO 
    … …  
    498498         WRITE(numout,*) '          Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) ) 
    499499 
    500          DO jj = 1, jpj 
    501             DO ji = 1, jpi 
    502                DO jk = mikt(ji,jj)+1, jpkm1 
    503                   zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s 
     500         DO jk = 2, jpkm1 
     501            DO jj = 1, jpj 
     502               DO ji = 1, jpi 
     503                  zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk)  !kz max = 300 cm2/s 
    504504               END DO 
    505505            END DO 
    … …  
    510510            DO jj = 1, jpj 
    511511               DO ji = 1, jpi 
    512                   ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 
     512                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
    513513               END DO 
    514514            END DO 
    … …  
    519519         DO jk = 1, jpk 
    520520            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zav_tide(:,:,jk)     * tmask_i(:,:) )   & 
    521                &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) ) 
     521               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 
    522522            ztpc = 1.E50 
    523523            DO jj = 1, jpj 
    … …  
    540540            END DO 
    541541            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   & 
    542                &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) ) 
     542               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 
    543543            WRITE(numout,*) '                jk= ', jk,'   ', ze_z * 1.e4,' cm2/s' 
    544544         END DO 
    … …  
    546546            zkz(:,:) = az_tmx(:,:,jk) /rn_n2min 
    547547            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   & 
    548                &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) ) 
     548               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 
    549549            WRITE(numout,*)  
    550550            WRITE(numout,*) '          N2 min - jk= ', jk,'   ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4,   & 

  • trunk/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r5118 r5120  
    342342         WRITE(numout,*) '                       NEMO team' 
    343343         WRITE(numout,*) '            Ocean General Circulation Model' 
    344          WRITE(numout,*) '                  version 3.4  (2011) ' 
     344         WRITE(numout,*) '                  version 3.6  (2015) ' 
    345345         WRITE(numout,*) 
    346346         WRITE(numout,*) 

  • trunk/NEMOGCM/NEMO/OPA_SRC/step.F90

    r5012 r5120  
    122122      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz 
    123123      IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value) 
    124          avt (:,:,:) = rn_avt0 * tmask(:,:,:) 
    125          avmu(:,:,:) = rn_avm0 * umask(:,:,:) 
    126          avmv(:,:,:) = rn_avm0 * vmask(:,:,:) 
     124         avt (:,:,:) = rn_avt0 * wmask (:,:,:) 
     125         avmu(:,:,:) = rn_avm0 * wumask(:,:,:) 
     126         avmv(:,:,:) = rn_avm0 * wvmask(:,:,:) 
    127127      ENDIF 
    128128      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths 
    … …  
    145145      ! 
    146146      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing 
    147                          CALL eos( tsb, rhd, gdept_0(:,:,:) )             ! before in situ density 
    148          IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    149             &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
    150             &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     147                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density 
     148         IF( ln_zps .AND. .NOT. ln_isfcav)                               & 
     149            &            CALL zps_hde    ( kstp, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
     150            &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     151         IF( ln_zps .AND.       ln_isfcav)                               & 
     152            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     153            &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     154            &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level 
    151155         IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator 
    152156                         CALL ldf_slp_grif( kstp ) 
    … …  
    177181          ! is necessary to compute momentum advection for the rhs of barotropic loop: 
    178182                            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 
    179             IF( ln_zps )    CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    180                &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
    181                &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     183            IF( ln_zps .AND. .NOT. ln_isfcav)                               & 
     184               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
     185               &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     186            IF( ln_zps .AND.       ln_isfcav)                               & 
     187               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     188               &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     189               &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    182190 
    183191                                  ua(:,:,:) = 0.e0             ! set dynamics trends to zero 
    … …  
    253261                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
    254262                             CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation 
    255          IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    256             &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
    257             &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     263            IF( ln_zps .AND. .NOT. ln_isfcav)                                & 
     264               &             CALL zps_hde    ( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
     265               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     266            IF( ln_zps .AND.       ln_isfcav)                                & 
     267               &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     268               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     269               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    258270      ELSE                                                  ! centered hpg  (eos then time stepping) 
    259271         IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case 
    260272                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation 
    261          IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    262          &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
    263          &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     273         IF( ln_zps .AND. .NOT. ln_isfcav)                                   & 
     274               &             CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
     275               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     276         IF( ln_zps .AND.       ln_isfcav)                                   &  
     277               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     278               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     279               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    264280         ENDIF 
    265281         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 

  • trunk/NEMOGCM/NEMO/OPA_SRC/timing.F90

    r3610 r5120  
    211211         WRITE(numtime,*) '                             NEMO team' 
    212212         WRITE(numtime,*) '                  Ocean General Circulation Model' 
    213          WRITE(numtime,*) '                        version 3.3  (2010) ' 
     213         WRITE(numtime,*) '                        version 3.6  (2015) ' 
    214214         WRITE(numtime,*) 
    215215         WRITE(numtime,*) '                        Timing Informations ' 
Note: See TracChangeset for help on using the changeset viewer.