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 10008 for NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90 – NEMO

Ignore:
Timestamp:
2018-07-26T17:58:13+02:00 (6 years ago)
Author:
acc
Message:

Branch: dev_r9956_ENHANCE05_ZAD_AIMP. Fix bug in trazdf.F90 to get closer to maintaining constant salinity in a simple GYRE test with uniform salinity. Also tidy dynzdf.F90 so results are the same with wi=0 and ln_zad_Aimp=T/F. See ticket #2042

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90

    r9979 r10008  
    7171      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars 
    7272      REAL(wp) ::   zzws, ze3va        !   -      - 
    73       REAL(wp) ::   z1_e3un, z1_e3vn   !   -      - 
     73      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      - 
    7474      REAL(wp) ::   zWu , zWv          !   -      - 
    7575      REAL(wp) ::   zWui, zWvi         !   -      - 
    … …  
    166166                  DO jj = 2, jpjm1  
    167167                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    168                         z1_e3un = 1._wp / e3u_n(ji,jj,jk) 
    169                         zzwi = ( ( avm   (ji+1,jj,jk  ) + avm   (ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
    170                            &     / e3uw_a(ji  ,jj,jk  ) ) * z1_e3un * wumask(ji,jj,jk  ) 
    171                         zzws = ( ( avm   (ji+1,jj,jk+1) + avm   (ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
    172                            &     / e3uw_a(ji  ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1) 
     168                        ze3ua   = e3u_a(ji,jj,jk) 
     169                        z1_e3ua = 1._wp / ze3ua 
     170                        zzwi = - zdt * ( avm   (ji+1,jj,jk  ) + avm   (ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     171                           &     / ( ze3ua * e3uw_n(ji  ,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     172                        zzws = - zdt * ( avm   (ji+1,jj,jk+1) + avm   (ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     173                           &     / ( ze3ua * e3uw_n(ji  ,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    173174                        zWu  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji+1,jj,jk  )   & 
    174175                           &               + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1)   ) 
    175                         zwi(ji,jj,jk) = - 1._wp * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un ) 
    176                         zws(ji,jj,jk) = - 1._wp * ( zzws - MAX( zWu, 0._wp ) * z1_e3un ) 
    177                         zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un ) 
     176                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWu, 0._wp ) * z1_e3ua 
     177                        zws(ji,jj,jk) = zzws + zdt * MAX( zWu, 0._wp ) * z1_e3ua 
     178                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3ua 
    178179                     END DO 
    179180                  END DO 
    … …  
    183184                  DO jj = 2, jpjm1  
    184185                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    185                         z1_e3un = 1._wp / e3u_n(ji,jj,jk) 
    186                         zzwi = ( ( avm   (ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   & 
    187                            &     / e3uw_a(ji  ,jj,jk  ) ) * z1_e3un * wumask(ji,jj,jk  ) 
    188                         zzws = ( ( avm   (ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   & 
    189                            &     / e3uw_a(ji  ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1) 
     186                        ze3ua   = e3u_a(ji,jj,jk) 
     187                        z1_e3ua = 1._wp / ze3ua 
     188                        zzwi = - zdt * ( avm   (ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   & 
     189                           &     / ( ze3ua * e3uw_n(ji  ,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     190                        zzws = - zdt * ( avm   (ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   & 
     191                           &     / ( ze3ua * e3uw_n(ji  ,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    190192                        zWu  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji+1,jj,jk  )   & 
    191193                           &               + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1)   ) 
    192                         zwi(ji,jj,jk) = - 1._wp * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un ) 
    193                         zws(ji,jj,jk) = - 1._wp * ( zzws - MAX( zWu, 0._wp ) * z1_e3un ) 
    194                         zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un ) 
     194                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWu, 0._wp ) * z1_e3ua  
     195                        zws(ji,jj,jk) = zzws + zdt * MAX( zWu, 0._wp ) * z1_e3ua  
     196                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3ua 
    195197                     END DO 
    196198                  END DO 
    … …  
    204206                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    205207                        ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
    206                         zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
    207                            &         / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    208                         zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
    209                            &         / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     208                        zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     209                           &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     210                        zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     211                           &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    210212                        zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) 
    211213                        zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 
    212                         zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWui, 0._wp ) )  
    213                         zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWus, 0._wp ) ) 
    214                         zwd(ji,jj,jk) = 1._wp + zzwi + zzws - zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     214                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWui, 0._wp )  
     215                        zws(ji,jj,jk) = zzws + zdt * MAX( zWus, 0._wp ) 
     216                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
    215217                     END DO 
    216218                  END DO 
    … …  
    221223                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    222224                        ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
    223                         zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    224                         zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     225                        zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     226                        zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    225227                        zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) 
    226228                        zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 
    227                         zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWui, 0._wp ) )  
    228                         zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWus, 0._wp ) ) 
    229                         zwd(ji,jj,jk) = 1._wp + zzwi + zzws - zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     229                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWui, 0._wp ) 
     230                        zws(ji,jj,jk) = zzws + zdt * MAX( zWus, 0._wp ) 
     231                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
    230232                     END DO 
    231233                  END DO 
    … …  
    240242                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    241243                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
    242                      zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
    243                         &         / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    244                      zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
    245                         &         / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    246                      zwi(ji,jj,jk) = - zdt * zzwi 
    247                      zws(ji,jj,jk) = - zdt * zzws 
    248                      zwd(ji,jj,jk) = 1._wp + zzwi + zzws 
     244                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     245                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     246                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     247                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     248                     zwi(ji,jj,jk) = zzwi 
     249                     zws(ji,jj,jk) = zzws 
     250                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    249251                  END DO 
    250252               END DO 
    … …  
    255257                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    256258                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
    257                      zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    258                      zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    259                      zwi(ji,jj,jk) = - zdt * zzwi 
    260                      zws(ji,jj,jk) = - zdt * zzws 
    261                      zwd(ji,jj,jk) = 1._wp + zzwi + zzws 
     259                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     260                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     261                     zwi(ji,jj,jk) = zzwi 
     262                     zws(ji,jj,jk) = zzws 
     263                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    262264                  END DO 
    263265               END DO 
    … …  
    361363                  DO jj = 2, jpjm1  
    362364                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    363                         z1_e3vn = 1._wp / e3v_n(ji,jj,jk) 
    364                         zzwi = ( ( avm   (ji,jj+1,jk  ) + avm   (ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
    365                            &     / e3vw_a(ji,jj  ,jk  ) ) * z1_e3vn * wvmask(ji,jj,jk  ) 
    366                         zzws = ( ( avm   (ji,jj+1,jk+1) + avm   (ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    367                            &     / e3vw_a(ji,jj  ,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1) 
     365                        ze3va   = e3v_a(ji,jj,jk) 
     366                        z1_e3va = 1._wp / ze3va 
     367                        zzwi = - zdt * ( avm   (ji,jj+1,jk  ) + avm   (ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     368                           &             / ( ze3va * e3vw_n(ji,jj  ,jk  ) ) * wvmask(ji,jj,jk  ) 
     369                        zzws = - zdt * ( avm   (ji,jj+1,jk+1) + avm   (ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     370                           &             / ( ze3va * e3vw_n(ji,jj  ,jk+1) ) * wvmask(ji,jj,jk+1) 
    368371                        zWv  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji,jj+1,jk  )   & 
    369372                           &               + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1)   ) 
    370                         zwi(ji,jj,jk) = - 1._wp * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn ) 
    371                         zws(ji,jj,jk) = - 1._wp * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn ) 
    372                         zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn ) 
     373                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWv, 0._wp ) * z1_e3va 
     374                        zws(ji,jj,jk) = zzws + zdt * MAX( zWv, 0._wp ) * z1_e3va 
     375                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3va 
    373376                     END DO 
    374377                  END DO 
    … …  
    378381                  DO jj = 2, jpjm1  
    379382                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    380                         z1_e3vn = 1._wp / e3v_n(ji,jj,jk) 
    381                         zzwi = ( ( avm   (ji,jj+1,jk  ) + avm(ji,jj,jk  ) )   & 
    382                            &     / e3vw_a(ji,jj  ,jk  ) ) * z1_e3vn * wvmask(ji,jj,jk  ) 
    383                         zzws = ( ( avm   (ji,jj+1,jk+1) + avm(ji,jj,jk+1) )   & 
    384                            &     / e3vw_a(ji  ,jj,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1) 
     383                        ze3va   = e3v_a(ji,jj,jk) 
     384                        z1_e3va = 1._wp / ze3va 
     385                        zzwi = - zdt * ( avm   (ji,jj+1,jk  ) + avm(ji,jj,jk  ) )   & 
     386                           &             / ( ze3va * e3vw_n(ji,jj  ,jk  ) ) * wvmask(ji,jj,jk  ) 
     387                        zzws = - zdt * ( avm   (ji,jj+1,jk+1) + avm(ji,jj,jk+1) )   & 
     388                           &             / ( ze3va * e3vw_n(ji  ,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    385389                        zWv  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji,jj+1,jk  )   & 
    386390                           &               + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1)   ) 
    387                         zwi(ji,jj,jk) = - 1._wp * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn ) 
    388                         zws(ji,jj,jk) = - 1._wp * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn ) 
    389                         zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn ) 
     391                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWv, 0._wp ) * z1_e3va 
     392                        zws(ji,jj,jk) = zzws + zdt * MAX( zWv, 0._wp ) * z1_e3va 
     393                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3va 
    390394                     END DO 
    391395                  END DO 
    … …  
    399403                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    400404                        ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at U-point 
    401                         zzwi = ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
    402                            &         / ( ze3va * e3vw_a(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    403                         zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    404                            &         / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    405                         zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) 
    406                         zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) 
    407                         zwi(ji,jj,jk) = - 1._wp * ( zzwi +  MIN( zWvi, 0._wp ) )  
    408                         zws(ji,jj,jk) = - 1._wp * ( zzws -  MAX( zWvs, 0._wp ) ) 
    409                         zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     405                        zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     406                           &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     407                        zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     408                           &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     409                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  ) 
     410                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) 
     411                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWvi, 0._wp ) 
     412                        zws(ji,jj,jk) = zzws + zdt * MAX( zWvs, 0._wp ) 
     413                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
    410414                     END DO 
    411415                  END DO 
    … …  
    416420                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    417421                        ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at U-point 
    418                         zzwi = ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_a(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    419                         zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    420                         zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) 
    421                         zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) 
    422                         zwi(ji,jj,jk) = - 1._wp * ( zzwi +  MIN( zWvi, 0._wp ) )  
    423                         zws(ji,jj,jk) = - 1._wp * ( zzws -  MAX( zWvs, 0._wp ) ) 
    424                         zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     422                        zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     423                        zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     424                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  ) 
     425                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) 
     426                        zwi(ji,jj,jk) = zzwi  - zdt * MIN( zWvi, 0._wp ) 
     427                        zws(ji,jj,jk) = zzws  + zdt * MAX( zWvs, 0._wp ) 
     428                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
    425429                     END DO 
    426430                  END DO 
    … …  
    439443                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    440444                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    441                      zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
    442                      zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
     445                     zwi(ji,jj,jk) = zzwi 
     446                     zws(ji,jj,jk) = zzws 
    443447                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    444448                  END DO 
    … …  
    452456                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    453457                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    454                      zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
    455                      zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
     458                     zwi(ji,jj,jk) = zzwi 
     459                     zws(ji,jj,jk) = zzws 
    456460                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    457461                  END DO 
Note: See TracChangeset for help on using the changeset viewer.