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 5189 for branches – NEMO

Changeset 5189 for branches


Ignore:
Timestamp:
2015-03-31T19:58:23+02:00 (9 years ago)
Author:
mathiot
Message:

ISF cleaning branch: simplification and bug correction in hpg_isf, zps_hde_isf, mixed layer definition, slope, diffusion, vertical advection and top friction

Location:
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM
Files:
20 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/CONFIG/ISOMIP/EXP00/iodef.xml

    r4924 r5189  
    3232      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE."/> <!-- 1d files --> 
    3333      <file_group id="3d" output_freq="3d"  output_level="10" enabled=".TRUE."/> <!-- 3d files -->     
    34       <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE.">  <!-- 5d files -->   
     34      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE."/>  <!-- 5d files -->   
    3535      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files --> 
    3636   <file id="file1" output_freq="1mo" name_suffix="_grid_T" description="ocean T grid variables" > 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90

    r5131 r5189  
    549549         &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    550550      IF( ln_zps .AND.        ln_isfcav)                            & 
    551          &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
    552          &                                        rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    553          &                                 gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level 
     551         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
     552         &                                        rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
    554553 
    555554      rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r5120 r5189  
    750750               &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level 
    751751            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 
     752               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,    &    ! Partial steps for top cell (ISF) 
     753               &                                  rhd, gru , grv , grui, grvi     )    ! of t, s, rd at the last ocean level 
    755754 
    756755#if defined key_zdfkpp 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r5123 r5189  
    646646      CALL iom_close( inum ) 
    647647       
    648 ! need to be define for the extended grid south of -80S 
    649 ! some point are undefined but you need to have e1 and e2 .NE. 0 
    650       WHERE (e1t==0.0_wp) 
    651          e1t=1.0e2 
    652       END WHERE 
    653       WHERE (e1v==0.0_wp) 
    654          e1v=1.0e2 
    655       END WHERE 
    656       WHERE (e1u==0.0_wp) 
    657          e1u=1.0e2 
    658       END WHERE 
    659       WHERE (e1f==0.0_wp) 
    660          e1f=1.0e2 
    661       END WHERE 
    662       WHERE (e2t==0.0_wp) 
    663          e2t=1.0e2 
    664       END WHERE 
    665       WHERE (e2v==0.0_wp) 
    666          e2v=1.0e2 
    667       END WHERE 
    668       WHERE (e2u==0.0_wp) 
    669          e2u=1.0e2 
    670       END WHERE 
    671       WHERE (e2f==0.0_wp) 
    672          e2f=1.0e2 
    673       END WHERE 
    674         
    675648    END SUBROUTINE hgr_read 
    676649     

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r5120 r5189  
    263263         END DO 
    264264      END DO 
    265       ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point 
     265      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point 
    266266      DO jj = 1, jpjm1 
    267267         DO ji = 1, fs_jpim1   ! vector loop 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5120 r5189  
    15381538  
    15391539 ! remove single point "bay" on isf coast line in the ice shelf draft' 
    1540          DO jk = 1, jpk 
     1540         DO jk = 2, jpk 
    15411541            WHERE (misfdep==0) misfdep=jpk 
    15421542            zmask=0 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5120 r5189  
    141141            &                                            rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    142142         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 
     143            &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,    &    ! Partial steps for top cell (ISF) 
     144            &                                            rhd, gru , grv , grui, grvi     )    ! of t, s, rd at the last ocean level 
    146145#endif 
    147146         !    

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5120 r5189  
    4444   USE wrk_nemo        ! Memory Allocation 
    4545   USE timing          ! Timing 
     46   USE iom 
    4647 
    4748   IMPLICIT NONE 
    … …  
    130131      INTEGER ::   ioptio = 0      ! temporary integer 
    131132      INTEGER ::   ios             ! Local integer output status for namelist read 
     133      !! 
     134      INTEGER  ::   ji, jj, jk, ikt    ! dummy loop indices      ISF 
     135      REAL(wp) ::   znad 
     136      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop, zrhd ! hypothesys on isf density 
     137      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_isf  ! density at bottom of ISF 
     138      REAL(wp), POINTER, DIMENSION(:,:)     ::  ziceload     ! density at bottom of ISF 
    132139      !! 
    133140      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     & 
    … …  
    191198      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 
    192199      !  
    193       ! initialisation of ice load 
    194       riceload(:,:)=0.0 
     200      ! initialisation of ice shelf load 
     201      IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0 
     202      IF (       ln_isfcav ) THEN 
     203         CALL wrk_alloc( jpi,jpj, 2,  ztstop)  
     204         CALL wrk_alloc( jpi,jpj,jpk, zrhd  ) 
     205         CALL wrk_alloc( jpi,jpj,     zrhdtop_isf, ziceload)  
     206         ! 
     207         IF(lwp) WRITE(numout,*) 
     208         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 
     209         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'   
     210 
     211         ! To use density and not density anomaly 
     212         znad=1._wp 
     213          
     214         ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 
     215         ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 
     216 
     217         ! compute density of the water displaced by the ice shelf  
     218         DO jk = 1, jpk 
     219            CALL eos(ztstop(:,:,:),fsdept_n(:,:,jk),zrhd(:,:,jk)) 
     220         END DO 
     221       
     222         ! compute rhd at the ice/oce interface (ice shelf side) 
     223         CALL eos(ztstop,risfdep,zrhdtop_isf) 
     224 
     225         ! Surface value + ice shelf gradient 
     226         ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 
     227         ziceload = 0._wp 
     228         DO jj = 1, jpj 
     229            DO ji = 1, jpi   ! vector opt. 
     230               ikt=mikt(ji,jj) 
     231               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 
     232               DO jk=2,ikt-1 
     233                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 
     234                     &                              * (1._wp - tmask(ji,jj,jk)) 
     235               END DO 
     236               IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 
     237                                  &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 
     238            END DO 
     239         END DO 
     240         riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5 
     241 
     242         CALL wrk_dealloc( jpi,jpj, 2,  ztstop)  
     243         CALL wrk_dealloc( jpi,jpj,jpk, zrhd  ) 
     244         CALL wrk_dealloc( jpi,jpj,     zrhdtop_isf, ziceload)  
     245      END IF 
    195246      ! 
    196247   END SUBROUTINE dyn_hpg_init 
    … …  
    465516      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    466517      !! 
    467       INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j                 ! dummy loop indices 
    468       REAL(wp) ::   zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1             ! temporary scalars 
    469       REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj, zrhd 
     518      INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j   ! dummy loop indices 
     519      REAL(wp) ::   zcoef0, zuap, zvap, znad                    ! temporary scalars 
     520      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj 
    470521      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop 
    471       REAL(wp), POINTER, DIMENSION(:,:)     ::  ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 
    472       !!---------------------------------------------------------------------- 
    473       ! 
    474       CALL wrk_alloc( jpi,jpj, 2, ztstop)  
    475       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 
    476       CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj)  
    477       ! 
    478      IF( kt == nit000 ) THEN 
    479          IF(lwp) WRITE(numout,*) 
    480          IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 
    481          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used' 
    482       ENDIF 
    483  
     522      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_oce, zpshpi, zpshpj 
     523      !!---------------------------------------------------------------------- 
     524      ! 
     525      CALL wrk_alloc( jpi,jpj,  2, ztstop)  
     526      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 
     527      CALL wrk_alloc( jpi,jpj,     zrhdtop_oce, zpshpi, zpshpj)  
     528      ! 
    484529      ! Local constant initialization 
    485530      zcoef0 = - grav * 0.5_wp 
     531    
    486532      ! To use density and not density anomaly 
    487 !      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume 
    488 !      ELSE                 ;     znad = 0._wp         ! Fixed volume 
    489 !      ENDIF 
    490533      znad=1._wp 
     534 
    491535      ! iniitialised to 0. zhpi zhpi  
    492536      zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 
    493537 
    494 !==================================================================================      
    495 !=====Compute iceload and contribution of the half first wet layer =================  
    496 !=================================================================================== 
    497  
    498       ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 
    499       ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 
    500  
    501       ! compute density of the water displaced by the ice shelf  
    502       zrhd = rhd ! save rhd 
    503       DO jk = 1, jpk 
    504            zdept(:,:)=gdept_1d(jk) 
    505            CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk)) 
    506       END DO 
    507       WHERE ( tmask(:,:,:) == 1._wp) 
    508         rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd 
    509       END WHERE 
    510        
    511       ! compute rhd at the ice/oce interface (ice shelf side) 
    512       CALL eos(ztstop,risfdep,zrhdtop_isf) 
    513  
    514       ! compute rhd at the ice/oce interface (ocean side) 
    515       DO ji=1,jpi 
    516         DO jj=1,jpj 
    517           ikt=mikt(ji,jj) 
    518           ztstop(ji,jj,1)=tsn(ji,jj,ikt,1) 
    519           ztstop(ji,jj,2)=tsn(ji,jj,ikt,2) 
    520         END DO 
    521       END DO 
    522       CALL eos(ztstop,risfdep,zrhdtop_oce) 
    523       ! 
    524       ! Surface value + ice shelf gradient 
    525       ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 
    526       ziceload = 0._wp 
    527       DO jj = 1, jpj 
    528          DO ji = 1, jpi   ! vector opt. 
    529             ikt=mikt(ji,jj) 
    530             ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 
    531             DO jk=2,ikt-1 
    532                ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 
    533                   &                              * (1._wp - tmask(ji,jj,jk)) 
    534             END DO 
    535             IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 
    536                                &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 
    537          END DO 
    538       END DO 
    539       riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5 
    540       ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 
     538      ! compute rhd at the ice/oce interface (ocean side) usefull to reduce unrealistic water current if homogene water 
     539!  TO BE DISCUSS do we compute density at w point for the ocean/ice atmosphere 
     540!      DO ji=1,jpi 
     541!        DO jj=1,jpj 
     542!          ikt=mikt(ji,jj) 
     543!          ztstop(ji,jj,1)=tsn(ji,jj,ikt,1) 
     544!          ztstop(ji,jj,2)=tsn(ji,jj,ikt,2) 
     545!        END DO 
     546!      END DO 
     547!      CALL eos(ztstop,risfdep,zrhdtop_oce) 
     548 
    541549      DO jj = 2, jpjm1 
    542550         DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    544552            ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 
    545553            ! we assume ISF is in isostatic equilibrium 
    546             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj  ,iktp1i)                                    & 
    547                &                                  * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   & 
    548                &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    & 
    549                &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
    550                &                                  + ( ziceload(ji+1,jj) - ziceload(ji,jj))                              )  
    551             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji  ,jj+1,iktp1j)                                    & 
    552                &                                  * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   & 
    553                &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    &  
    554                &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
    555                &                                  + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             )  
     554            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,iktp1i) * ( znad + rhd(ji+1,jj  ,iktp1i) )   & 
     555               &                                  - fse3w(ji  ,jj  ,ikt   ) * ( znad + rhd(ji  ,jj  ,ikt   ) )   & 
     556               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                       )  
     557            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,iktp1j) * ( znad + rhd(ji  ,jj+1,iktp1j) )   & 
     558               &                                  - fse3w(ji  ,jj  ,ikt   ) * ( znad + rhd(ji  ,jj  ,ikt   ) )   & 
     559               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj) )                      )  
     560!            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj  ,iktp1i)                                    & 
     561!               &                                  * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   & 
     562!               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    & 
     563!               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
     564!               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                              )  
     565!            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji  ,jj+1,iktp1j)                                    & 
     566!               &                                  * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   & 
     567!               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    &  
     568!               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
     569!               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj) )                             )  
    556570            ! s-coordinate pressure gradient correction (=0 if z coordinate) 
    557571            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    … …  
    564578         END DO 
    565579      END DO 
    566 !==================================================================================      
    567 !===== Compute partial cell contribution for the top cell =========================  
    568 !================================================================================== 
    569       DO jj = 2, jpjm1 
    570          DO ji = fs_2, fs_jpim1   ! vector opt. 
    571             iku = miku(ji,jj) ;  
    572             zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp 
    573             ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
    574             ! u direction 
    575             IF ( iku .GT. 1 ) THEN 
    576                ! case iku 
    577                zhpi(ji,jj,iku)   =  zcoef0 / e1u(ji,jj) * ze3wu                                            & 
    578                       &                                 * ( rhd    (ji+1,jj,iku) + rhd   (ji,jj,iku)       & 
    579                       &                                   + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad ) 
    580                ! corrective term ( = 0 if z coordinate ) 
    581                zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj) 
    582                ! zhpi will be added in interior loop 
    583                ua(ji,jj,iku)     = ua(ji,jj,iku) + zuap 
    584                ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure   
    585                IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj)  = zhpi(ji,jj,iku) 
    586  
    587                ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps) 
    588                zhpiint        =  zcoef0 / e1u(ji,jj)                                                               &     
    589                   &           * (  fse3w(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          & 
    590                   &                                         + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   & 
    591                   &              - fse3w(ji  ,jj  ,iku+1) * ( (rhd(ji  ,jj,iku+1) + znad)                          & 
    592                   &                                         + (rhd(ji  ,jj,iku  ) + znad) ) * tmask(ji  ,jj,iku)   ) 
    593                zhpi(ji,jj,iku+1) =  zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint  
    594             END IF 
    595                 
    596             ! v direction 
    597             ikv = mikv(ji,jj) 
    598             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)) 
    599             IF ( ikv .GT. 1 ) THEN 
    600                ! case ikv 
    601                zhpj(ji,jj,ikv)   =  zcoef0 / e2v(ji,jj) * ze3wv                                            & 
    602                      &                                  * ( rhd(ji,jj+1,ikv) + rhd   (ji,jj,ikv)           & 
    603                      &                                    + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 
    604                ! corrective term ( = 0 if z coordinate ) 
    605                zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj) 
    606                ! zhpi will be added in interior loop 
    607                va(ji,jj,ikv)      = va(ji,jj,ikv) + zvap 
    608                ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure   
    609                IF (mbkv(ji,jj) == ikv + 1)  zpshpj(ji,jj)  =  zhpj(ji,jj,ikv)  
    610                 
    611                ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps) 
    612                zhpjint        =  zcoef0 / e2v(ji,jj)                                                              & 
    613                   &           * (  fse3w(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         & 
    614                   &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)    & 
    615                   &              - fse3w(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         & 
    616                   &                                       + (rhd(ji,jj  ,ikv  ) + znad) ) * tmask(ji,jj  ,ikv)  ) 
    617                zhpj(ji,jj,ikv+1) =  zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 
    618             END IF 
    619          END DO 
    620       END DO 
    621580 
    622581!==================================================================================      
    623582!===== Compute interior value =====================================================  
    624583!================================================================================== 
    625  
    626       DO jj = 2, jpjm1 
    627          DO ji = fs_2, fs_jpim1   ! vector opt. 
    628             iku=miku(ji,jj); ikv=mikv(ji,jj) 
    629             DO jk = 2, jpkm1 
     584      ! interior value (2=<jk=<jpkm1) 
     585      DO jk = 2, jpkm1 
     586         DO jj = 2, jpjm1 
     587            DO ji = fs_2, fs_jpim1   ! vector opt. 
    630588               ! hydrostatic pressure gradient along s-surfaces 
    631                ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc) 
    632                zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                              & 
    633                   &                            + zcoef0 / e1u(ji,jj)                                                           & 
    634                   &                                     * ( fse3w(ji+1,jj  ,jk) * ( (rhd(ji+1,jj,jk  ) + znad)                 & 
    635                   &                                                     + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1)   & 
    636                   &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                 & 
    637                   &                                                     + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   )  
     589               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
     590                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   & 
     591                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   ) 
     592               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
     593                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   & 
     594                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   ) 
    638595               ! s-coordinate pressure gradient correction 
    639                ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)  
    640                zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                    & 
    641                   &            * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 
    642                ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
    643  
    644                ! hydrostatic pressure gradient along s-surfaces 
    645                ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc) 
    646                zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                              & 
    647                   &                            + zcoef0 / e2v(ji,jj)                                                           & 
    648                   &                                     * ( fse3w(ji  ,jj+1,jk) * ( (rhd(ji,jj+1,jk  ) + znad)                 & 
    649                   &                                                     + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1)   & 
    650                   &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                 & 
    651                   &                                                     + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   ) 
    652                ! s-coordinate pressure gradient correction 
    653                ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
    654                zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                     & 
    655                   &            * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 
     596               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     597                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
     598               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     599                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
    656600               ! add to the general momentum trend 
    657                va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 
     601               ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
     602               va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 
    658603            END DO 
    659604         END DO 
    660605      END DO 
    661  
    662 !==================================================================================      
    663 !===== Compute bottom cell contribution (partial cell) ============================  
    664 !================================================================================== 
    665  
    666       DO jj = 2, jpjm1 
    667          DO ji = 2, jpim1 
    668             iku = mbku(ji,jj) 
    669             ikv = mbkv(ji,jj) 
    670  
    671             IF (iku .GT. 1) THEN 
    672                ! remove old value (interior case) 
    673                zuap            = -zcoef0 * ( rhd   (ji+1,jj  ,iku) + rhd   (ji,jj,iku) + 2._wp * znad )   & 
    674                      &                   * ( fsde3w(ji+1,jj  ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 
    675                ua(ji,jj,iku)   = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 
    676                ! put new value 
    677                ! -zpshpi to avoid double contribution of the partial step in the top layer  
    678                zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  / e1u(ji,jj) 
    679                zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj)  
    680                ua(ji,jj,iku)   =  ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 
    681             END IF 
    682             ! v direction 
    683             IF (ikv .GT. 1) THEN 
    684                ! remove old value (interior case) 
    685                zvap            = -zcoef0 * ( rhd   (ji  ,jj+1,ikv) + rhd   (ji,jj,ikv) + 2._wp * znad )   & 
    686                      &                   * ( fsde3w(ji  ,jj+1,ikv) - fsde3w(ji,jj,ikv) )   / e2v(ji,jj) 
    687                va(ji,jj,ikv)   = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 
    688                ! put new value 
    689                ! -zpshpj to avoid double contribution of the partial step in the top layer  
    690                zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj)     / e2v(ji,jj) 
    691                zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)    
    692                va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 
    693             END IF 
    694          END DO 
    695       END DO 
    696       
    697       ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 
    698       rhd = zrhd 
    699       ! 
    700       CALL wrk_dealloc( jpi,jpj,2, ztstop) 
    701       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 
    702       CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 
     606      ! 
     607      CALL wrk_dealloc( jpi,jpj,2  , ztstop) 
     608      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj) 
     609      CALL wrk_dealloc( jpi,jpj    , zrhdtop_oce, zpshpi, zpshpj) 
    703610      ! 
    704611   END SUBROUTINE hpg_isf 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r5120 r5189  
    8282      ENDIF 
    8383       
    84       DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical 
     84      !                             ! Vertical momentum advection at uw and vw-points 
     85      DO jk = 2, jpk                 
    8586         DO jj = 2, jpj                   ! vertical fluxes  
    86             DO ji = fs_2, jpi             ! vector opt. 
    87                zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 
     87            DO ji = fs_2, jpi 
     88               zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
    8889            END DO 
    8990         END DO 
    90          DO jj = 2, jpjm1                 ! vertical momentum advection at w-point 
    91             DO ji = fs_2, fs_jpim1        ! vector opt. 
    92                zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) ) 
    93                zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) ) 
     91         DO jj = 2, jpjm1                 ! interior vertical momentum advection at w-point 
     92            DO ji = fs_2, fs_jpim1 
     93               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) 
     94               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) 
    9495            END DO   
    9596         END DO    
    9697      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 
     98      DO jj = 2, jpjm1                    ! First/last level advective fluxes 
     99         DO ji = fs_2, fs_jpim1           ! w is zero at k=1 and k=jpk 
     100            zwuw(ji,jj, 1 ) = 0._wp 
     101            zwvw(ji,jj, 1 ) = 0._wp 
     102            zwuw(ji,jj,jpk) = 0._wp 
     103            zwvw(ji,jj,jpk) = 0._wp 
     104         END DO 
     105      END DO 
    118106 
    119107      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points 
    … …  
    205193         DO jj = 2, jpj                    
    206194            DO ji = fs_2, jpi             ! vector opt. 
    207                zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 
     195               zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
    208196            END DO 
    209197         END DO 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r5120 r5189  
    162162            !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0 
    163163            !                                !    else  tmask(ik+1)  = 0   =>   pn2(ik+1) = 0   =>   zdzr divides by 1 
    164             !                                !          umask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2 
     164            !                                !          tmask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2 
    165165            !                                ! NB: 1/(tmask+1) = (1-.5*tmask)  substitute a / by a *  ==> faster 
    166166            zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp )              & 
    … …  
    188188            DO jj = 2, jpjm1 
    189189               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 
    197             ENDDO 
     190                  zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp) - 0.5 * ( risfdep(ji,jj) + risfdep(ji+1,jj  ) ) 
     191                  zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp) - 0.5 * ( risfdep(ji,jj) + risfdep(ji  ,jj+1) ) 
     192               END DO 
     193            END DO 
    198194         ELSE 
    199195            DO jj = 2, jpjm1 
    … …  
    201197                  zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp) 
    202198                  zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp) 
    203                ENDDO 
    204             ENDDO 
     199               END DO 
     200            END DO 
    205201         END IF 
     202! PM + GM can be optimized by using : zuslp_hmlpu(ji,jj)= uslpml(ji,jj) / zhmlpu (ji,jj) 
     203 
    206204         DO jk = 2, jpkm1                            !* Slopes at u and v points 
    207205            DO jj = 2, jpjm1 
    … …  
    220218                  zfj = MAX( omlmask(ji,jj,jk), omlmask(ji  ,jj+1,jk) ) 
    221219                  ! thickness of water column between surface and level k at u/v point 
    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)) ) 
     220                  zdepu = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji+1,jj  ,jk) )                            & 
     221                                   - ( risfdep(ji,jj)    + risfdep(ji+1,jj)    ) - fse3u(ji,jj,miku(ji,jj)) ) 
     222                  zdepv = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji,jj+1,jk) )                              & 
     223                                   - ( risfdep(ji,jj)    + risfdep(ji,jj+1)    ) - fse3v(ji,jj,mikv(ji,jj)) ) 
    226224                  ! 
    227225                  zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                          & 
    … …  
    315313                  !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 
    316314                  zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0 
    317                   zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp ) 
     315                  zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - fsdepw(ji,jj,mikt(ji,jj)), 10._wp ) 
    318316                  zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) & 
    319317                     &            + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk) 
    … …  
    426424         DO jj = 2, jpjm1  
    427425            DO ji = fs_2, fs_jpim1   ! vector opt.  
    428                uslp(ji,jj,1) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,1) - fsdept_b(ji ,jj ,1) )  * umask(ji,jj,1)  
    429                vslp(ji,jj,1) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,1) - fsdept_b(ji ,jj ,1) )  * vmask(ji,jj,1)  
     426               uslp(ji,jj,1)  = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,1) - fsdept_b(ji ,jj ,1) )  * umask(ji,jj,1)  
     427               vslp(ji,jj,1)  = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,1) - fsdept_b(ji ,jj ,1) )  * vmask(ji,jj,1)  
    430428               wslpi(ji,jj,1) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,1) - fsdepw_b(ji-1,jj,1) ) * tmask(ji,jj,1) * 0.5  
    431429               wslpj(ji,jj,1) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,1) - fsdepw_b(ji,jj-1,1) ) * tmask(ji,jj,1) * 0.5  
    … …  
    436434            DO jj = 2, jpjm1  
    437435               DO ji = fs_2, fs_jpim1   ! vector opt.  
    438                   uslp(ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)  
    439                   vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)  
     436                  uslp(ji,jj,jk)  = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)  
     437                  vslp(ji,jj,jk)  = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)  
    440438                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) & 
    441439                    &                              * wmask(ji,jj,jk) * 0.5  
    … …  
    752750            DO ji = 1, jpi 
    753751               ik = nmln(ji,jj) - 1 
    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 
     752               IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 
     753               ELSE                ; omlmask(ji,jj,jk) = 0._wp 
    758754               ENDIF 
    759755            END DO 
    … …  
    777773            ! 
    778774            !                        !- vertical density gradient for u- and v-slopes (from dzr at T-point) 
    779             iku = MIN(  MAX( miku(ji,jj)+1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1) 
    780             ikv = MIN(  MAX( mikv(ji,jj)+1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   ! 
     775            iku = MIN(  MAX( nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1) 
     776            ikv = MIN(  MAX( nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   ! 
    781777            zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) ) 
    782778            zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) ) 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r5147 r5189  
    851851               pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
    852852                  &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
    853                   &            / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
     853                  &            / fse3w(ji,jj,jk) * wmask(ji,jj,jk) 
    854854            END DO 
    855855         END DO 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r5149 r5189  
    108108      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
    109109      REAL(wp) ::  zcoef0, zbtr, ztra            !   -      - 
    110       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d 
    111       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdkt, zdk1t, zdit, zdjt, ztfw  
     110      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zdkt, zdk1t, z2d 
     111      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdit, zdjt , ztfw  
    112112      !!---------------------------------------------------------------------- 
    113113      ! 
    114114      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso') 
    115115      ! 
    116       CALL wrk_alloc( jpi, jpj,      z2d )  
    117       CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t )  
     116      CALL wrk_alloc( jpi, jpj,      zdkt, zdk1t, z2d )  
     117      CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt , ztfw )  
    118118      ! 
    119119 
    … …  
    168168         !!   II - horizontal trend  (full) 
    169169         !!---------------------------------------------------------------------- 
    170 !!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t )  
     170!CDIR PARALLEL DO PRIVATE( zdk1t )  
     171         !                                                ! =============== 
     172         DO jk = 1, jpkm1                                 ! Horizontal slab 
     173            !                                             ! =============== 
    171174            ! 1. Vertical tracer gradient at level jk and jk+1 
    172175            ! ------------------------------------------------ 
    173          !  
    174          ! interior value  
    175          DO jk = 2, jpkm1                
    176             DO jj = 1, jpj 
    177                DO ji = 1, jpi   ! vector opt. 
    178                   zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn  ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) 
    179                   ! 
    180                   zdkt(ji,jj,jk)  = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn  ) ) * wmask(ji,jj,jk) 
    181                END DO 
    182             END DO 
    183          END DO 
    184          ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
    185          zdk1t(:,:,1) = ( ptb(:,:,1,jn  ) - ptb(:,:,2,jn) ) * wmask(:,:,2) 
    186          zdkt (:,:,1) = zdk1t(:,:,1) 
    187          IF ( ln_isfcav ) THEN 
    188             DO jj = 1, jpj 
    189                DO ji = 1, jpi   ! vector opt. 
    190                   ikt = mikt(ji,jj) ! surface level 
    191                   zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 
    192                   zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 
    193                END DO 
    194             END DO 
    195          END IF 
    196  
    197          ! 2. Horizontal fluxes 
    198          ! --------------------    
    199          DO jk = 1, jpkm1 
     176            ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
     177            zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1) 
     178            ! 
     179            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:) 
     180            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 
     181            ENDIF 
     182 
     183            ! 2. Horizontal fluxes 
     184            ! --------------------    
    200185            DO jj = 1 , jpjm1 
    201186               DO ji = 1, fs_jpim1   ! vector opt. 
    … …  
    203188                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
    204189                  ! 
    205                   zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   & 
    206                      &             + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1. ) 
    207                   ! 
    208                   zmskv = 1. / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   & 
    209                      &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. ) 
     190                  zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
     191                     &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     192                  ! 
     193                  zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   & 
     194                     &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. ) 
    210195                  ! 
    211196                  zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
    … …  
    213198                  ! 
    214199                  zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
    215                      &              + zcof1 * (  zdkt (ji+1,jj,jk) + zdk1t(ji,jj,jk)      & 
    216                      &                         + zdk1t(ji+1,jj,jk) + zdkt (ji,jj,jk)  )  ) * umask(ji,jj,jk) 
     200                     &              + zcof1 * (  zdkt (ji+1,jj  ) + zdk1t(ji,jj)      & 
     201                     &                         + zdk1t(ji+1,jj  ) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
    217202                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
    218                      &              + zcof2 * (  zdkt (ji,jj+1,jk) + zdk1t(ji,jj,jk)      & 
    219                      &                         + zdk1t(ji,jj+1,jk) + zdkt (ji,jj,jk)  )  ) * vmask(ji,jj,jk)                   
     203                     &              + zcof2 * (  zdkt (ji  ,jj+1) + zdk1t(ji,jj)      & 
     204                     &                         + zdk1t(ji  ,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                   
    220205               END DO 
    221206            END DO 
    … …  
    322307      END DO 
    323308      ! 
    324       CALL wrk_dealloc( jpi, jpj, z2d )  
    325       CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t )  
     309      CALL wrk_dealloc( jpi, jpj,      zdkt, zdk1t, z2d )  
     310      CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt , ztfw )  
    326311      ! 
    327312      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso') 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r5120 r5189  
    222222      ! 
    223223      IF( nn_isf > 0 ) THEN 
    224          zfact = 0.5e0 
     224         zfact = 0.5_wp 
    225225         DO jj = 2, jpj 
    226226            DO ji = fs_2, fs_jpim1 
    … …  
    235235               ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    236236!                  zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 
    237                   zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 
     237                  zt_frz = -1.9_wp !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 
    238238               ! compute trend 
    239                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          & 
    240                      &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          & 
    241                      &               - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) & 
     239                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                & 
     240                     &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)               & 
     241                     &                     - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) & 
    242242                     &           * r1_hisf_tbl(ji,jj) 
    243                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          & 
    244                      &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) 
     243                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                                & 
     244                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal))               & 
     245                     &           * r1_hisf_tbl(ji,jj) 
    245246               END DO 
    246247    
    … …  
    248249               ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    249250!               zpress = grav*rau0*fsdept(ji,jj,ikb)*1.e-04 
    250                zt_frz = -1.9 !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress ) 
     251               zt_frz = -1.9_wp !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress ) 
    251252               ! compute trend 
    252                tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           & 
    253                   &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          & 
    254                   &                  - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) &  
     253               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 & 
     254                  &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)               & 
     255                  &                        - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) &  
    255256                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
    256                tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           & 
    257                   &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)  
     257               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                                 & 
     258                  &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal))               & 
     259                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)  
     260 
    258261            END DO 
    259262         END DO 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r5120 r5189  
    196196   END SUBROUTINE zps_hde 
    197197   ! 
    198    SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv,   & 
    199       &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  & 
    200       &                   pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 
     198   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  & 
     199      &                          prd, pgru, pgrv, pgrui, pgrvi ) 
    201200      !!---------------------------------------------------------------------- 
    202201      !!                     ***  ROUTINE zps_hde  *** 
    … …  
    252251      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    253252      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom) 
    254       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru, pmrv      ! hor. sum  of prd at u- & v-pts (bottom) 
    255       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom) 
    256       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom) 
    257253      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) 
    261254      ! 
    262255      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    … …  
    269262      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf') 
    270263      ! 
    271       pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
     264      pgtu (:,:,:)=0.0_wp ; pgtv(:,:,:) =0.0_wp ; 
    272265      pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ; 
    273       zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
    274       zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
     266      zti  (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
     267      zhi  (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
    275268      ! 
    276269      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
    … …  
    280273               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    281274               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
     275               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku) 
     276               ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv) 
     277               ! 
     278               ! i- direction 
     279               IF( ze3wu >= 0._wp ) THEN      ! case 1 
     280                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
     281                  ! interpolated values of tracers 
     282                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     283                  ! gradient of  tracers 
     284                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     285               ELSE                           ! case 2 
     286                  zmaxu = -ze3wu / fse3w(ji,jj,iku) 
     287                  ! interpolated values of tracers 
     288                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     289                  ! gradient of tracers 
     290                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     291               ENDIF 
     292               ! 
     293               ! j- direction 
     294               IF( ze3wv >= 0._wp ) THEN      ! case 1 
     295                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
     296                  ! interpolated values of tracers 
     297                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     298                  ! gradient of tracers 
     299                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     300               ELSE                           ! case 2 
     301                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
     302                  ! interpolated values of tracers 
     303                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     304                  ! gradient of tracers 
     305                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     306               ENDIF 
     307            END DO 
     308         END DO 
     309         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     310         ! 
     311      END DO 
     312 
     313      ! horizontal derivative of density anomalies (rd) 
     314      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
     315         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
     316         ! 
     317         DO jj = 1, jpjm1 
     318            DO ji = 1, jpim1 
     319               iku = mbku(ji,jj) 
     320               ikv = mbkv(ji,jj) 
     321               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku) 
     322               ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv) 
     323 
     324               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
     325               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
     326               ENDIF 
     327               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1 
     328               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2 
     329               ENDIF 
     330 
     331            END DO 
     332         END DO 
     333 
     334         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
     335         ! step and store it in  zri, zrj for each  case 
     336         CALL eos( zti, zhi, zri ) 
     337         CALL eos( ztj, zhj, zrj ) 
     338 
     339         ! Gradient of density at the last level  
     340         DO jj = 1, jpjm1 
     341            DO ji = 1, jpim1 
     342               iku = mbku(ji,jj) 
     343               ikv = mbkv(ji,jj) 
     344               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku) 
     345               ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv) 
     346 
     347               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
     348               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2 
     349               ENDIF 
     350               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1 
     351               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2 
     352               ENDIF 
     353            END DO 
     354         END DO 
     355         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
     356         ! 
     357      END IF 
     358         ! (ISH)  compute grui and gruvi 
     359      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
     360         DO jj = 1, jpjm1 
     361            DO ji = 1, jpim1 
     362               iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1 
     363               ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1 
     364               ! 
    282365               ! (ISF) case partial step top and bottom in adjacent cell in vertical 
    283366               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 
    284367               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
    285368               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
    286                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    287                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    288                ! 
     369               ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku) 
     370               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)  
    289371               ! i- direction 
    290372               IF( ze3wu >= 0._wp ) THEN      ! case 1 
    291                   zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
    292                   ! interpolated values of tracers 
    293                   zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     373                  zmaxu = ze3wu / fse3w(ji+1,jj,iku+1) 
     374                  ! interpolated values of tracers 
     375                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 
     376                  ! gradient of tracers 
     377                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     378               ELSE                           ! case 2 
     379                  zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 
     380                  ! interpolated values of tracers 
     381                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 
    294382                  ! gradient of  tracers 
    295                   pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    296                ELSE                           ! case 2 
    297                   zmaxu = -ze3wu / fse3w(ji,jj,iku) 
    298                   ! interpolated values of tracers 
    299                   zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
    300                   ! gradient of tracers 
    301                   pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     383                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    302384               ENDIF 
    303385               ! 
    304386               ! j- direction 
    305387               IF( ze3wv >= 0._wp ) THEN      ! case 1 
    306                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
    307                   ! interpolated values of tracers 
    308                   ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
    309                   ! gradient of tracers 
    310                   pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    311                ELSE                           ! case 2 
    312                   zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
    313                   ! interpolated values of tracers 
    314                   ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
    315                   ! gradient of tracers 
    316                   pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    317                ENDIF 
    318             END DO 
    319          END DO 
    320          CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     388                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1) 
     389                  ! interpolated values of tracers 
     390                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 
     391                  ! gradient of tracers 
     392                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     393               ELSE                           ! case 2 
     394                  zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1) 
     395                  ! interpolated values of tracers 
     396                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 
     397                  ! gradient of tracers 
     398                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     399               ENDIF 
     400            END DO!! 
     401         END DO!! 
     402         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    321403         ! 
    322404      END DO 
    … …  
    324406      ! horizontal derivative of density anomalies (rd) 
    325407      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    326          pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
    327          pgzu(:,:)=0.0_wp   ; pgzv(:,:)=0.0_wp ; 
    328          pmru(:,:)=0.0_wp   ; pmru(:,:)=0.0_wp ; 
    329          pge3ru(:,:)=0.0_wp ; pge3rv(:,:)=0.0_wp ; 
    330          DO jj = 1, jpjm1 
    331             DO ji = 1, jpim1 
    332                iku = mbku(ji,jj) 
    333                ikv = mbkv(ji,jj) 
    334                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    335                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    336  
    337                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1 
    338                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2 
    339                ENDIF 
    340                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv    ! j-direction: case 1 
    341                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) + ze3wv    ! -     -      case 2 
    342                ENDIF 
    343             END DO 
    344          END DO 
    345           
     408         pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ; 
     409         DO jj = 1, jpjm1 
     410            DO ji = 1, jpim1 
     411               iku = miku(ji,jj) 
     412               ikv = mikv(ji,jj) 
     413               ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku) 
     414               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)  
     415 
     416               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
     417               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
     418               ENDIF 
     419 
     420               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1 
     421               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2 
     422               ENDIF 
     423 
     424            END DO 
     425         END DO 
     426 
    346427         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    347428         ! step and store it in  zri, zrj for each  case 
    … …  
    352433         DO jj = 1, jpjm1 
    353434            DO ji = 1, jpim1 
    354                iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    355                ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    356                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    357                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    358                IF( ze3wu >= 0._wp ) THEN  
    359                   pgzu(ji,jj) = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku) 
    360                   pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
    361                   pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 1  
    362                   pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  & 
    363                                 * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji  ,jj    ) + prd(ji+1,jj,ikum1) + 2._wp) & 
    364                                    - fse3w(ji  ,jj,iku)          * ( prd(ji  ,jj,iku) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2 
    365                ELSE   
    366                   pgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu) 
    367                   pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2 
    368                   pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) )   ! i: 2 
    369                   pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  & 
    370                                 * (  fse3w(ji+1,jj,iku)          * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) & 
    371                                    -(fse3w(ji  ,jj,iku) + ze3wu) * ( zri(ji  ,jj    ) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2 
    372                ENDIF 
    373                IF( ze3wv >= 0._wp ) THEN 
    374                   pgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv)  
    375                   pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1 
    376                   pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )   ! j: 1 
    377                   pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  & 
    378                                 * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj      ) + prd(ji,jj+1,ikvm1) + 2._wp) & 
    379                                    - fse3w(ji,jj  ,ikv)          * ( prd(ji,jj  ,ikv) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2 
    380                ELSE  
    381                   pgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv) 
    382                   pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2 
    383                   pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )   ! j: 2 
    384                   pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  & 
    385                                 * (  fse3w(ji,jj+1,ikv)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) & 
    386                                    -(fse3w(ji,jj  ,ikv) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2 
    387                ENDIF 
    388             END DO 
    389          END DO 
    390          CALL lbc_lnk( pgru   , 'U', -1. )   ;   CALL lbc_lnk( pgrv   , 'V', -1. )   ! Lateral boundary conditions 
    391          CALL lbc_lnk( pmru   , 'U',  1. )   ;   CALL lbc_lnk( pmrv   , 'V',  1. )   ! Lateral boundary conditions 
    392          CALL lbc_lnk( pgzu   , 'U', -1. )   ;   CALL lbc_lnk( pgzv   , 'V', -1. )   ! Lateral boundary conditions 
    393          CALL lbc_lnk( pge3ru , 'U', -1. )   ;   CALL lbc_lnk( pge3rv , 'V', -1. )   ! Lateral boundary conditions 
    394          ! 
    395       END IF 
    396          ! (ISH)  compute grui and gruvi 
    397       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
    398          DO jj = 1, jpjm1 
    399             DO ji = 1, jpim1 
    400                iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1 
    401                ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1 
    402                ! 
    403                ! (ISF) case partial step top and bottom in adjacent cell in vertical 
    404                ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 
    405                ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
    406                ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
    407                ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))  
    408                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)) 
    409                ! i- direction 
    410                IF( ze3wu >= 0._wp ) THEN      ! case 1 
    411                   zmaxu = ze3wu / fse3w(ji+1,jj,iku+1) 
    412                   ! interpolated values of tracers 
    413                   zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 
    414                   ! gradient of tracers 
    415                   pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    416                ELSE                           ! case 2 
    417                   zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 
    418                   ! interpolated values of tracers 
    419                   zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 
    420                   ! gradient of  tracers 
    421                   pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    422                ENDIF 
    423                ! 
    424                ! j- direction 
    425                IF( ze3wv >= 0._wp ) THEN      ! case 1 
    426                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1) 
    427                   ! interpolated values of tracers 
    428                   ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 
    429                   ! gradient of tracers 
    430                   pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    431                ELSE                           ! case 2 
    432                   zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1) 
    433                   ! interpolated values of tracers 
    434                   ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 
    435                   ! gradient of tracers 
    436                   pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    437                ENDIF 
    438             END DO!! 
    439          END DO!! 
    440          CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    441          ! 
    442       END DO 
    443  
    444       ! horizontal derivative of density anomalies (rd) 
    445       IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    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 ; 
    450  
    451          DO jj = 1, jpjm1 
    452             DO ji = 1, jpim1 
    453                iku = miku(ji,jj) 
    454                ikv = mikv(ji,jj) 
    455                ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
    456                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)) 
    457  
    458                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
    459                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
    460                ENDIF 
    461                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1 
    462                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2 
    463                ENDIF 
    464             END DO 
    465          END DO 
    466  
    467          ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    468          ! step and store it in  zri, zrj for each  case 
    469          CALL eos( zti, zhi, zri )   
    470          CALL eos( ztj, zhj, zrj ) 
    471  
    472          ! Gradient of density at the last level  
    473          DO jj = 1, jpjm1 
    474             DO ji = 1, jpim1 
    475435               iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 
    476436               ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 
    477                ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
    478                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)) 
    479                IF( ze3wu >= 0._wp ) THEN 
    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)                                                                  & 
    484                                 * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   & 
    485                                    - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1 
    486                ELSE 
    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)                                                                   & 
    491                                 * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  & 
    492                                    -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2 
    493                ENDIF 
    494                IF( ze3wv >= 0._wp ) THEN 
    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)                                                                  &  
    499                                 * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  & 
    500                                    - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1 
    501                                   ! + 2 due to the formulation in density and not in anomalie in hpg sco 
    502                ELSE 
    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)                                                                   & 
    507                                 * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 
    508                                    -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2 
    509                ENDIF 
     437               ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku) 
     438               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)  
     439 
     440               IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) )          ! i: 1 
     441               ELSE                      ; pgrui(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj  ,iku) - zri(ji,jj    ) )      ! i: 2 
     442               ENDIF 
     443 
     444               IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji  ,jj      ) - prd(ji,jj,ikv) ) ! j: 1 
     445               ELSE                      ; pgrvi(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji  ,jj+1,ikv) - zrj(ji,jj    ) ) ! j: 2 
     446               ENDIF 
     447 
    510448            END DO 
    511449         END DO 
    512450         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 
    516451         ! 
    517452      END IF   

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r5120 r5189  
    171171            END DO 
    172172         END DO 
     173         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition 
    173174         IF ( ln_isfcav ) THEN 
    174175            DO jj = 2, jpjm1 
    … …  
    202203               END DO 
    203204            END DO 
     205            CALL lbc_lnk( tfrua, 'U', 1. )   ;   CALL lbc_lnk( tfrva, 'V', 1. )      ! Lateral boundary condition 
    204206         END IF 
    205207         ! 
    206          CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition 
    207208         ! 
    208209         IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        & 
    … …  
    295296         bfrua(:,:) = - bfrcoef2d(:,:) 
    296297         bfrva(:,:) = - bfrcoef2d(:,:) 
     298         ! 
     299         IF ( ln_isfcav ) THEN 
     300            IF(ln_tfr2d) THEN 
     301               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
     302               CALL iom_open('tfr_coef.nc',inum) 
     303               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array 
     304               CALL iom_close(inum) 
     305               tfrcoef2d(:,:) = rn_tfri1 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 
     306            ELSE 
     307               tfrcoef2d(:,:) = rn_tfri1  ! initialize tfrcoef2d to the namelist variable 
     308            ENDIF 
     309            ! 
     310            tfrua(:,:) = - tfrcoef2d(:,:) 
     311            tfrva(:,:) = - tfrcoef2d(:,:) 
     312         END IF 
    297313         ! 
    298314      CASE( 2 ) 
    … …  
    336352            bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable 
    337353         ENDIF 
     354          
     355         IF ( ln_isfcav ) THEN 
     356            IF(ln_tfr2d) THEN 
     357               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
     358               CALL iom_open('tfr_coef.nc',inum) 
     359               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! bfrcoef2d is used as tmp array 
     360               CALL iom_close(inum) 
     361               ! 
     362               tfrcoef2d(:,:) = rn_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 
     363            ELSE 
     364               tfrcoef2d(:,:) = rn_tfri2  ! initialize tfrcoef2d to the namelist variable 
     365            ENDIF 
     366         END IF 
    338367         ! 
    339368         IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all 
    … …  
    346375               END DO 
    347376            END DO 
     377            IF ( ln_isfcav ) THEN 
     378               DO jj = 1, jpj 
     379                  DO ji = 1, jpi 
     380                     ikbt = mikt(ji,jj) 
     381                     ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_tfrz0 ))**2._wp 
     382                     tfrcoef2d(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 
     383                     tfrcoef2d(ji,jj) = MIN(tfrcoef2d(ji,jj), rn_tfri2_max) 
     384                  END DO 
     385               END DO 
     386            END IF 
    348387         ENDIF 
    349388         ! 
    … …  
    398437             zminbfr = MIN(  zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) )  ) 
    399438             zmaxbfr = MAX(  zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) )  ) 
     439! (ISF) 
     440             IF ( ln_isfcav ) THEN 
     441                ikbu = miku(ji,jj)       ! deepest ocean level at u- and v-points 
     442                ikbv = mikv(ji,jj) 
     443                zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt 
     444                zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt 
     445                IF( ABS( tfrcoef2d(ji,jj) ) > zfru ) THEN 
     446                   IF( ln_ctl ) THEN 
     447                      WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu 
     448                      WRITE(numout,*) 'BFR ', ABS( tfrcoef2d(ji,jj) ), zfru 
     449                   ENDIF 
     450                   ictu = ictu + 1 
     451                ENDIF 
     452                IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN 
     453                   IF( ln_ctl ) THEN 
     454                      WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbv 
     455                      WRITE(numout,*) 'BFR ', tfrcoef2d(ji,jj), zfrv 
     456                   ENDIF 
     457                   ictv = ictv + 1 
     458                ENDIF 
     459                zmintfr = MIN(  zmintfr, MIN( zfru, ABS( tfrcoef2d(ji,jj) ) )  ) 
     460                zmaxtfr = MAX(  zmaxtfr, MIN( zfrv, ABS( tfrcoef2d(ji,jj) ) )  ) 
     461             END IF 
     462! END ISF 
    400463         END DO 
    401464      END DO 
    … …  
    405468         CALL mpp_min( zminbfr ) 
    406469         CALL mpp_max( zmaxbfr ) 
     470         IF ( ln_isfcav) CALL mpp_min( zmintfr ) 
     471         IF ( ln_isfcav) CALL mpp_max( zmaxtfr ) 
    407472      ENDIF 
    408473      IF( .NOT.ln_bfrimp) THEN 
    … …  
    411476         WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points ' 
    412477         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr 
    413          WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr 
    414478         WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary' 
     479         IF ( ln_isfcav ) WRITE(numout,*) ' Top friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr 
    415480      ENDIF 
    416481      ENDIF 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r4990 r5189  
    3131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m] 
    3232   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    33    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
     33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: depth of the last T-point inside the mixed layer 
    3434 
    3535   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    … …  
    111111      END DO 
    112112      ! 
    113       ! w-level of the turbocline 
     113      ! w-level of the turbocline and mixing layer (iom_use) 
    114114      imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point 
    115115      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10  
    116116         DO jj = 1, jpj 
    117117            DO ji = 1, jpi 
    118                imkt = mikt(ji,jj) 
    119                IF( avt (ji,jj,jk) < avt_c )   imld(ji,jj) = MAX( imkt, jk )      ! Turbocline  
     118               IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline  
    120119            END DO 
    121120         END DO 
    … …  
    127126            iikn = nmln(ji,jj) 
    128127            imkt = mikt(ji,jj) 
    129             hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )  * ssmask(ji,jj)    ! Turbocline depth  
    130             hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth 
    131             hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt )            )  * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer 
     128            hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth  
     129            hmlp (ji,jj) = fsdepw(ji,jj,iikn ) * ssmask(ji,jj)    ! Mixed layer depth 
     130            hmlpt(ji,jj) = fsdept(ji,jj,iikn-1) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer 
    132131         END DO 
    133132      END DO 
    134133      IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode 
    135          CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth 
    136          CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth 
     134         IF ( iom_use("mldr10_1") ) THEN 
     135            IF( .NOT. ln_isfcav ) CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth 
     136            IF(       ln_isfcav ) CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness 
     137         END IF 
     138         IF ( iom_use("mldkz5") ) THEN 
     139            IF( .NOT. ln_isfcav ) CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth 
     140            IF(       ln_isfcav ) CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness 
     141         END IF 
    137142      ENDIF 
    138143       

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r4990 r5189  
    4545   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gtsu, gtsv   !: horizontal gradient of T, S bottom u-point  
    4646   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gru , grv    !: horizontal gradient of rd at bottom u-point 
    47    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aru , arv     
    48    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gzu , gzv     
    49    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ge3ru, ge3rv   !: horizontal gradient of T, S and rd at top v-point   
    5047 
    5148   !! (ISF) interpolated gradient (only used for ice shelf case)  
    … …  
    5350   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gtui, gtvi   !: horizontal gradient of T, S and rd at top u-point  
    5451   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   grui, grvi   !: horizontal gradient of T, S and rd at top v-point   
    55    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   arui, arvi   !: horizontal average  of rd          at top v-point   
    56    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gzui, gzvi   !: horizontal gradient of z           at top v-point   
    57    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ge3rui, ge3rvi   !: horizontal gradient of T, S and rd at top v-point   
    5852   !! (ISF) ice load 
    5953   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   riceload 
    … …  
    10498         &     spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       & 
    10599         &     gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     & 
    106          &     aru(jpi,jpj)      , arv(jpi,jpj)      ,                     & 
    107          &     gzu(jpi,jpj)      , gzv(jpi,jpj)      ,                     & 
    108100         &     gru(jpi,jpj)      , grv(jpi,jpj)      ,                     & 
    109          &     ge3ru(jpi,jpj)    , ge3rv(jpi,jpj)    ,                     & 
    110101         &     gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts),                     & 
    111          &     arui(jpi,jpj)     , arvi(jpi,jpj)     ,                     & 
    112          &     gzui(jpi,jpj)     , gzvi(jpi,jpj)     ,                     & 
    113          &     ge3rui(jpi,jpj)   , ge3rvi(jpi,jpj)   ,                     & 
    114102         &     grui(jpi,jpj)     , grvi(jpi,jpj)     ,                     & 
    115103         &     riceload(jpi,jpj),                             STAT=ierr(2) ) 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/step.F90

    r5147 r5189  
    150150            &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    151151         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 
     152            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
     153            &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
    155154         IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator 
    156155                         CALL ldf_slp_grif( kstp ) 
    … …  
    181180          ! is necessary to compute momentum advection for the rhs of barotropic loop: 
    182181                            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 
    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 
     182            IF( ln_zps .AND. .NOT. ln_isfcav )                               & 
     183               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient 
     184               &                                          rhd, gru , grv     )  ! of t, s, rd at the last ocean level 
     185            IF( ln_zps .AND.       ln_isfcav )                                          & 
     186               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
     187               &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
    190188 
    191189                                  ua(:,:,:) = 0.e0             ! set dynamics trends to zero 
    … …  
    266264               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    267265            IF( ln_zps .AND.       ln_isfcav)                                & 
    268                &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
    269                &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    270                &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
     266               &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
     267               &                                           rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
    271268      ELSE                                                  ! centered hpg  (eos then time stepping) 
    272269         IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case 
    … …  
    276273               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    277274         IF( ln_zps .AND.       ln_isfcav)                                   &  
    278                &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
    279                &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    280                &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
     275               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
     276               &                                           rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
    281277         ENDIF 
    282278         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90

    r5120 r5189  
    8686            &            CALL zps_hde    ( kstp, jptra, trn, gtru, gtrv )   ! Partial steps: now horizontal gradient of passive 
    8787         IF( ln_zps .AND.        ln_isfcav)        & 
    88             &            CALL zps_hde_isf( kstp, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! Partial steps: now horizontal gradient of passive 
     88            &            CALL zps_hde_isf( kstp, jptra, trn, gtru, gtrv, gtrui, gtrvi )  ! Partial steps: isf case 
    8989                                                                ! tracers at the bottom ocean level 
    9090         ! 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/TOP_SRC/trcini.F90

    r5120 r5189  
    146146        &    CALL zps_hde    ( nit000, jptra, trn, gtru, gtrv  )  ! Partial steps: before horizontal gradient 
    147147      IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav )   & 
    148         &    CALL zps_hde_isf( nit000, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )       ! tracers at the bottom ocean level 
     148        &    CALL zps_hde_isf( nit000, jptra, trn, gtru, gtrv, gtrui, gtrvi )       ! tracers at the bottom ocean level 
    149149 
    150150 
Note: See TracChangeset for help on using the changeset viewer.