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

Changeset 1380


Ignore:
Timestamp:
2009-04-06T13:00:29+02:00 (15 years ago)
Author:
rblod
Message:

Change scale factor computation in VVL case, see ticket #397

Location:
branches/dev_004_VVL/NEMO/OPA_SRC
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domvvl.F90

    r1362 r1380  
    1414   !!   dom_vvl         : defined scale factors & depths at each time step 
    1515   !!   dom_vvl_ini     : defined coefficients to distribute ssh on each layers 
    16    !!   dom_vvl_ssh     : defined the ocean sea level at each time step 
    1716   !!---------------------------------------------------------------------- 
    1817   !! * Modules used 
    … …  
    3231   !! * Routine accessibility 
    3332   PUBLIC dom_vvl_ini    ! called by dom_init.F90 
    34    PUBLIC dom_vvl_ssh    ! called by trazdf.F90 
    3533   PUBLIC dom_vvl        ! called by istate.F90 and step.F90 
    36    PUBLIC dom_vvl_sf_ini !  
    37    PUBLIC dom_vvl_sf     !  
    3834 
    3935   !! * Module variables 
     36   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hu_0, hv_0 
     37   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ee_t, ee_u, ee_v, ee_f   !: 
     38 
    4039   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !: 
    4140      mut, muu, muv, muf                            !: 
    … …  
    6564      !! 
    6665      !!---------------------------------------------------------------------- 
    67       INTEGER  ::   ji, jj, jk, zpk 
    68       REAL(wp), DIMENSION(jpi,jpj) ::   zmut, zmuu, zmuv, zmuf   ! 2D workspace 
     66      INTEGER  ::   ji, jj, jk 
     67      REAL(wp) ::   zcoefu, zcoefv, zcoeff 
    6968      !!---------------------------------------------------------------------- 
    7069 
    … …  
    7473         WRITE(numout,*) '~~~~~~~~~~~   compute coef. used to spread ssh over each layers' 
    7574      ENDIF 
    76  
    77       IF( ln_zps )  CALL ctl_stop( 'dom_vvl_ini : option ln_zps is incompatible with variable volume option key_vvl') 
    7875 
    7976#if defined key_zco  ||  defined key_dynspg_rl 
    … …  
    9390 
    9491      ! mu computation 
    95       zmut(:,:)   = 0.e0 
    96       zmuu(:,:)   = 0.e0 
    97       zmuv(:,:)   = 0.e0 
    98       zmuf(:,:)   = 0.e0 
    99       mut (:,:,:) = 0.e0 
    100       muu (:,:,:) = 0.e0 
    101       muv (:,:,:) = 0.e0 
    102       muf (:,:,:) = 0.e0 
    103  
    104       DO jj = 1, jpj 
    105          DO ji = 1, jpi 
    106             zpk = mbathy(ji,jj) - 1 
    107             DO jk = 1, zpk 
    108                zmut(ji,jj) = zmut(ji,jj) + fse3t_0(ji,jj,jk) * SUM( fse3t_0(ji,jj,jk:zpk) ) 
    109                zmuu(ji,jj) = zmuu(ji,jj) + fse3u_0(ji,jj,jk) * SUM( fse3u_0(ji,jj,jk:zpk) ) 
    110                zmuv(ji,jj) = zmuv(ji,jj) + fse3v_0(ji,jj,jk) * SUM( fse3v_0(ji,jj,jk:zpk) ) 
    111                zmuf(ji,jj) = zmuf(ji,jj) + fse3f_0(ji,jj,jk) * SUM( fse3f_0(ji,jj,jk:zpk) ) 
    112             END DO 
    113          END DO 
    114       END DO 
    115  
    116       DO jj = 1, jpj 
    117          DO ji = 1, jpi 
    118             zpk = mbathy(ji,jj) - 1 
    119             DO jk = 1, zpk 
    120 #if defined key_sigma_vvl 
    121                mut(ji,jj,jk) = 1./SUM( fse3t_0(ji,jj,1:zpk) )  
    122                muu(ji,jj,jk) = 1./SUM( fse3u_0(ji,jj,1:zpk) )  
    123                muv(ji,jj,jk) = 1./SUM( fse3v_0(ji,jj,1:zpk) )  
    124                muf(ji,jj,jk) = 1./SUM( fse3f_0(ji,jj,1:zpk) )  
    125 #else 
    126                mut(ji,jj,jk) = SUM( fse3t_0(ji,jj,jk:zpk) ) / zmut(ji,jj) 
    127                muu(ji,jj,jk) = SUM( fse3u_0(ji,jj,jk:zpk) ) / zmuu(ji,jj) 
    128                muv(ji,jj,jk) = SUM( fse3v_0(ji,jj,jk:zpk) ) / zmuv(ji,jj) 
    129                muf(ji,jj,jk) = SUM( fse3f_0(ji,jj,jk:zpk) ) / zmuf(ji,jj) 
    130 #endif 
    131             END DO 
    132          END DO 
    133       END DO 
     92      ! -------------- 
     93      ! define ee_t, u, v and f as in sigma coordinate (ee_t = 1/ht, ...) 
     94      ee_t(:,:) = fse3t_0(:,:,1)        ! Lower bound : thickness of the first model level 
     95      ee_u(:,:) = fse3u_0(:,:,1) 
     96      ee_v(:,:) = fse3v_0(:,:,1) 
     97      ee_f(:,:) = fse3f_0(:,:,1) 
     98      DO jk = 2, jpkm1                   ! Sum of the masked vertical scale factors 
     99         ee_t(:,:) = ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 
     100         ee_u(:,:) = ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 
     101         ee_v(:,:) = ee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 
     102         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask 
     103            ee_f(:,jj) = ee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk) 
     104         END DO 
     105      END DO   
     106      !                                  ! Compute and mask the inverse of the local depth at T, U, V and F points 
     107      ee_t(:,:) = 1. / ee_t(:,:) * tmask(:,:,1) 
     108      ee_u(:,:) = 1. / ee_u(:,:) * umask(:,:,1) 
     109      ee_v(:,:) = 1. / ee_v(:,:) * vmask(:,:,1) 
     110      DO jj = 1, jpjm1                         ! f-point case fmask cannot be used  
     111         ee_f(:,jj) = 1. / ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 
     112      END DO 
     113      CALL lbc_lnk( ee_f, 'F', 1. )       ! lateral boundary condition on ee_f 
     114      ! 
     115      DO jk = 1, jpk 
     116         mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk)   ! at T levels 
     117         muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk)   ! at T levels 
     118         muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk)   ! at T levels 
     119      END DO 
     120      DO jk = 1, jpk 
     121         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask 
     122               muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels 
     123         END DO 
     124         muf(:,jpj,jk) = 0.e0 
     125      END DO 
     126      CALL lbc_lnk( muf, 'F', 1. )       ! lateral boundary condition on ee_f 
     127 
     128 
     129!!debug print 
     130!         ii=50   ;   ij = 50 
     131!         do jk= 1, jpk  
     132!         WRITE(numout,*) 'domvvl GM  : h0=', SUM( fse3t_0(ii,ij,1:jk) * tmask(ii,ij,1:jk) ), 'e3t0=', fse3t_0(ii,ij,jk),   & 
     133!            &            'e3t =', fse3t_0(ii,ij,jk) * ( 1 +  mut(ii,ij,jk) ), 'mut', mut(ii,ij,jk),   & 
     134!            &            'h =', SUM( fse3t_0(ii,ij,1:jk) * ( 1 +  mut(ii,ij,1:jk) ) * tmask(ii,ij,1:jk) ) 
     135!         end do 
     136!!end debug print 
     137 
     138      ! Reference ocean depth at U- and V-points 
     139      hu_0(:,:) = 0.e0     
     140      hv_0(:,:) = 0.e0 
     141      DO jk = 1, jpk 
     142         hu_0(:,:) = hu_0(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 
     143         hv_0(:,:) = hv_0(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 
     144      END DO 
     145 
     146      ! before and now Sea Surface Height at u-, v-, f-points 
     147      DO jj = 1, jpjm1 
     148         DO ji = 1, jpim1 
     149            zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
     150            zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
     151            zcoeff = 0.25 * fmask(ji,jj,1) 
     152!!gm bug used of fmask, even if thereafter multiplied by muf  which is correctly masked) 
     153            sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
     154               &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )    
     155            sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
     156               &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )    
     157            sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 & 
     158               &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )                
     159            sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
     160               &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )    
     161            sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &  
     162               &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )      
     163            sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 & 
     164               &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )                
     165         END DO 
     166      END DO 
     167      ! Boundaries conditions 
     168      CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
     169      CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
     170      CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. ) 
     171      ! 
    134172 
    135173 
    136174   END SUBROUTINE dom_vvl_ini 
    137  
    138175 
    139176 
    … …  
    147184      !! * Local declarations 
    148185      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    149       REAL(wp), DIMENSION(jpi,jpj) :: zsshf    ! 2D workspace 
    150       !!---------------------------------------------------------------------- 
    151  
    152       ! Sea Surface Height at u-v-fpoints 
    153       DO jj = 1, jpjm1 
    154          DO ji = 1,jpim1 
    155             sshu(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )  & 
    156                &          * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
    157                &            + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
    158             ! 
    159             sshv(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )  & 
    160                &          * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
    161                &            + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
    162             ! 
    163             zsshf(ji,jj) = 0.25 * fmask(ji,jj,1)                 & 
    164                &           * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   & 
    165                &             + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) 
    166          END DO 
    167       END DO 
    168  
    169       ! Boundaries conditions 
    170       CALL lbc_lnk( sshu,  'U', 1. ) 
    171       CALL lbc_lnk( sshv,  'V', 1. ) 
    172       CALL lbc_lnk( zsshf, 'F', 1. ) 
     186      !!---------------------------------------------------------------------- 
     187 
     188 !     IF( kt == nit000 ) THEN 
     189 !        IF(lwp) WRITE(numout,*) 
     190 !        IF(lwp) WRITE(numout,*) 'dom_vvl : ' 
     191 !        IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
     192 !     ENDIF 
    173193 
    174194      ! Scale factors at T levels 
    175195      DO jk = 1, jpkm1 
    176          fse3t(:,:,jk) = fse3t_0(:,:,jk) * ( 1 + sshn(:,:)  * mut(:,:,jk) ) 
    177          fse3u(:,:,jk) = fse3u_0(:,:,jk) * ( 1 + sshu(:,:)  * muu(:,:,jk) ) 
    178          fse3v(:,:,jk) = fse3v_0(:,:,jk) * ( 1 + sshv(:,:)  * muv(:,:,jk) ) 
    179          fse3f(:,:,jk) = fse3f_0(:,:,jk) * ( 1 + zsshf(:,:) * muf(:,:,jk) ) 
    180       END DO 
    181  
    182       ! Scale factors at W levels 
    183       fse3w (:,:,1) = fse3t(:,:,1) 
    184       fse3uw(:,:,1) = fse3u(:,:,1) 
    185       fse3vw(:,:,1) = fse3v(:,:,1) 
    186       DO jk = 2, jpk 
    187          fse3w (:,:,jk) = 0.5 * ( fse3t(:,:,jk-1) + fse3t(:,:,jk) ) 
    188          fse3uw(:,:,jk) = 0.5 * ( fse3u(:,:,jk-1) + fse3u(:,:,jk) ) 
    189          fse3vw(:,:,jk) = 0.5 * ( fse3v(:,:,jk-1) + fse3v(:,:,jk) ) 
    190       END DO 
    191  
    192       ! T and W points depth 
    193       fsdept(:,:,1) = 0.5 * fse3w(:,:,1) 
    194       fsdepw(:,:,1) = 0.e0 
    195       fsde3w(:,:,1) = fsdept(:,:,1) - sshn(:,:) 
    196       DO jk = 2, jpk 
    197          fsdept(:,:,jk) = fsdept(:,:,jk-1) + fse3w(:,:,jk) 
    198          fsdepw(:,:,jk) = fsdepw(:,:,jk-1) + fse3t(:,:,jk-1) 
    199          fsde3w(:,:,jk) = fsdept(:,:,jk  ) - sshn(:,:) 
    200       END DO 
    201  
    202       IF( MINVAL(fse3t(:,:,:)) < 0.0 ) THEN 
    203          CALL ctl_stop('E R R O R : dom_vvl','Level thickness fse3t less than zero.') 
    204          nstop = nstop+1 
    205       ENDIF 
     196         fse3t (:,:,jk) = fse3t_n (:,:,jk) 
     197         fse3u (:,:,jk) = fse3u_n (:,:,jk) 
     198         fse3v (:,:,jk) = fse3v_n (:,:,jk) 
     199         fse3f (:,:,jk) = fse3f_n (:,:,jk) 
     200         fse3w (:,:,jk) = fse3w_n (:,:,jk) 
     201         fse3uw(:,:,jk) = fse3uw_n(:,:,jk) 
     202         fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 
     203      END DO 
     204 
     205      DO jk = 1, jpkm1 
     206         DO jj = 1, jpj 
     207            DO ji = 1, jpi 
     208               fsdept(ji,jj,jk) = fsdept_n(ji,jj,jk) 
     209               fsdepw(ji,jj,jk) = fsdepw_n(ji,jj,jk) 
     210               fsde3w(ji,jj,jk) = fsde3w_n(ji,jj,jk) 
     211            END DO 
     212         END DO 
     213      END DO 
    206214 
    207215      ! Local depth or Inverse of the local depth of the water column at u- and v-points 
    208216      ! ------------------------------ 
    209  
    210       ! Ocean depth at U- and V-points 
    211       hu(:,:) = 0.e0    ;    hv(:,:) = 0.e0 
    212  
    213       DO jk = 1, jpk 
    214          hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 
    215          hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 
    216       END DO 
    217  
    218  
    219       ! Inverse of the local depth 
    220       hur(:,:) = fse3u(:,:,1)             ! Lower bound : thickness of the first model level 
    221       hvr(:,:) = fse3v(:,:,1) 
    222        
    223       DO jk = 2, jpk                      ! Sum of the vertical scale factors 
    224          hur(:,:) = hur(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 
    225          hvr(:,:) = hvr(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 
    226       END DO 
    227  
    228       ! Compute and mask the inverse of the local depth 
    229       hur(:,:) = 1. / hur(:,:) * umask(:,:,1) 
    230       hvr(:,:) = 1. / hvr(:,:) * vmask(:,:,1) 
     217      hu(:,:) = hu_0(:,:) + sshu_n(:,:) 
     218      hv(:,:) = hv_0(:,:) + sshv_n(:,:) 
     219 
     220      ! masked  inverse of the local depth 
     221      hur(:,:) = 1. / MAX( hu(:,:), fse3u_0(:,:,1) ) * umask(:,:,1) 
     222      hvr(:,:) = 1. / MAX( hv(:,:), fse3v_0(:,:,1) ) * vmask(:,:,1) 
    231223 
    232224   END SUBROUTINE dom_vvl 
    233225 
    234  
    235  
    236    SUBROUTINE dom_vvl_ssh( kt )  
    237       !!---------------------------------------------------------------------- 
    238       !!                ***  ROUTINE dom_vvl_ssh  *** 
    239       !!                    
    240       !! ** Purpose :  compute the ssha field used in tra_zdf, dynnxt, tranxt  
    241       !!               and dynspg routines 
    242       !!---------------------------------------------------------------------- 
    243       !! * Arguments 
    244       INTEGER, INTENT( in )  ::    kt                         ! time step 
    245  
    246       !! * Local declarations 
    247       INTEGER  ::   ji, jj, jk                                ! dummy loop indices 
    248       REAL(wp) ::   z2dt, zraur                               ! temporary scalars 
    249       REAL(wp), DIMENSION(jpi,jpj) ::   zun, zvn, zhdiv       ! 2D workspace 
    250       !!---------------------------------------------------------------------- 
    251  
    252       !! Sea surface height after (ssha) in variable volume case 
    253       !! --------------------------====-----===============----- 
    254       !! ssha is needed for the tracer time-stepping (trazdf_imp or 
    255       !! tra_nxt), for the ssh time-stepping (dynspg_*) and for the 
    256       !! dynamics time-stepping (dynspg_flt or dyn_nxt, and wzv). 
    257       !! un, vn (or un_b and vn_b) and emp are needed, so it must be 
    258       !! done after the call of oce_sbc. 
    259  
    260       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000 
    261          r2dt(:) =  rdttra(:)                       ! = rdtra (restarting with Euler time stepping) 
    262       ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1 
    263          r2dt(:) = 2. * rdttra(:)                   ! = 2 rdttra (leapfrog) 
    264       ENDIF 
    265  
    266       z2dt  = r2dt(1) 
    267       zraur = 1. / rauw 
    268  
    269       ! Vertically integrated quantities 
    270       ! -------------------------------- 
    271       IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN 
    272          zun(:,:) = 0.e0    ;    zvn(:,:) = 0.e0 
    273          ! 
    274          DO jk = 1, jpkm1 
    275             !                                               ! Vertically integrated transports (now) 
    276             zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
    277             zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
    278          END DO 
    279       ELSE ! lk_dynspg_ts is true 
    280          zun (:,:) = un_b(:,:)    ;    zvn (:,:) = vn_b(:,:) 
    281       ENDIF 
    282  
    283       ! Horizontal divergence of barotropic transports 
    284       !-------------------------------------------------- 
    285       DO jj = 2, jpjm1 
    286          DO ji = 2, jpim1   ! vector opt. 
    287             zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )    & 
    288                &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )    & 
    289                &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )    & 
    290                &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )  & 
    291                &           / ( e1t(ji,jj) * e2t(ji,jj) ) 
    292          END DO 
    293       END DO 
    294  
    295 #if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts ) 
    296       ! open boundaries (div must be zero behind the open boundary) 
    297       !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    298       IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east 
    299       IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west 
    300       IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north 
    301       IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south 
    302 #endif 
    303  
    304 #if defined key_bdy 
    305       zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) 
    306 #endif 
    307  
    308       CALL lbc_lnk( zhdiv, 'T', 1. ) 
    309  
    310       ! Sea surface elevation time stepping 
    311       ! ----------------------------------- 
    312       IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN 
    313          ssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1) 
    314       ELSE ! lk_dynspg_ts is true 
    315          ssha(:,:) = (  sshb_b(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    316       ENDIF 
    317  
    318  
    319    END SUBROUTINE dom_vvl_ssh 
    320  
    321    SUBROUTINE  dom_vvl_sf( zssh, gridp, sfe3 ) 
    322       !!---------------------------------------------------------------------- 
    323       !!                ***  ROUTINE dom_vvl_sf  *** 
    324       !!                    
    325       !! ** Purpose :  compute vertical scale factor at each time step 
    326       !!---------------------------------------------------------------------- 
    327       !! * Arguments 
    328       CHARACTER(LEN=1)                , INTENT( in ) :: gridp   ! grid point type 
    329       REAL(wp), DIMENSION(jpi,jpj)    , INTENT( in ) :: zssh    ! 2D workspace 
    330       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: sfe3    ! 3D workspace 
    331  
    332       !! * Local declarations 
    333       INTEGER  ::   jk                                      ! dummy loop indices 
    334       !!---------------------------------------------------------------------- 
    335  
    336       SELECT CASE (gridp) 
    337  
    338       CASE ('U') 
    339          ! 
    340          DO jk = 1, jpk 
    341             sfe3(:,:,jk)  = fve3u_0(:,:,jk) * ( 1 + zssh(:,:) * muu(:,:,jk) ) 
    342          ENDDO 
    343  
    344       CASE ('V') 
    345          ! 
    346          DO jk = 1, jpk 
    347             sfe3(:,:,jk)  = fse3v_0(:,:,jk) * ( 1 + zssh(:,:) * muv(:,:,jk) ) 
    348          ENDDO 
    349  
    350       CASE ('T') 
    351          ! 
    352          DO jk = 1, jpk 
    353             sfe3(:,:,jk)  = fse3t_0(:,:,jk) * ( 1 + zssh(:,:) * mut(:,:,jk) ) 
    354          ENDDO 
    355  
    356       END SELECT 
    357  
    358    END SUBROUTINE dom_vvl_sf 
    359  
    360    SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini ) 
    361       !!---------------------------------------------------------------------- 
    362       !!                ***  ROUTINE dom_vvl_sf_ini  *** 
    363       !!                    
    364       !! ** Purpose :  affect the appropriate vertical scale factor. It is done 
    365       !!               to avoid compiling error when using key_zco cpp key 
    366       !!---------------------------------------------------------------------- 
    367       !! * Arguments 
    368       CHARACTER(LEN=1)                , INTENT( in ) :: gridp      ! grid point type 
    369       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini    ! 3D workspace 
    370       !!---------------------------------------------------------------------- 
    371  
    372       SELECT CASE (gridp) 
    373  
    374       CASE ('U') 
    375          ! 
    376          sfe3ini(:,:,:)  = fse3u(:,:,:) 
    377  
    378       CASE ('V') 
    379          ! 
    380          sfe3ini(:,:,:)  = fse3v(:,:,:) 
    381  
    382       CASE ('T') 
    383          ! 
    384          sfe3ini(:,:,:)  = fse3t(:,:,:) 
    385  
    386       END SELECT 
    387  
    388    END SUBROUTINE dom_vvl_sf_ini 
    389226#else 
    390227   !!---------------------------------------------------------------------- 
    … …  
    395232   SUBROUTINE dom_vvl 
    396233   END SUBROUTINE dom_vvl 
    397    SUBROUTINE dom_vvl_ssh( kt )  
    398       !! * Arguments 
    399       INTEGER, INTENT( in )  ::    kt        ! time step 
    400       WRITE(*,*) 'dom_vvl_ssh: You should not have seen this print! error?', kt 
    401    END SUBROUTINE dom_vvl_ssh 
    402    SUBROUTINE dom_vvl_sf( zssh, gridp, sfe3 )  
    403       !! * Arguments 
    404       CHARACTER(LEN=1)                , INTENT( in ) :: gridp   ! grid point type 
    405       REAL(wp), DIMENSION(jpi,jpj)    , INTENT( in ) :: zssh    ! 2D workspace 
    406       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3    ! 3D workspace 
    407       sfe3(:,:,:) = 0.e0 
    408       WRITE(*,*) 'sfe3: You should not have seen this print! error?', gridp 
    409       WRITE(*,*) 'sfe3: You should not have seen this print! error?', zssh(1,1) 
    410    END SUBROUTINE dom_vvl_sf 
    411    SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini )  
    412       !! * Arguments 
    413       CHARACTER(LEN=1)                , INTENT( in ) :: gridp    ! grid point type 
    414       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini  ! 3D workspace 
    415       sfe3ini(:,:,:) = 0.e0 
    416       WRITE(*,*) 'sfe3ini: You should not have seen this print! error?', gridp 
    417    END SUBROUTINE dom_vvl_sf_ini 
    418234#endif 
    419235 

  • branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r1368 r1380  
    188188         CALL lbc_lnk( zsshua, 'U', 1. )    ;    CALL lbc_lnk( zsshva, 'V', 1. ) 
    189189 
    190          ! Scale factors at before and after time step 
    191          ! ------------------------------------------- 
    192          CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ;    CALL dom_vvl_sf( zsshua, 'U', zfse3ua ) 
    193          CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ;    CALL dom_vvl_sf( zsshva, 'V', zfse3va ) 
    194  
    195  
    196190         ! Thickness weighting 
    197191         ! ------------------- 

  • branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r1368 r1380  
    529529            ! Boundaries conditions 
    530530            CALL lbc_lnk( zsshun_e, 'U', 1. )    ;    CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
    531  
    532             ! Scale factors at before and after time step 
    533             ! ------------------------------------------- 
    534             CALL dom_vvl_sf( zsshun_e, 'U', zfse3un_e ) ;   CALL dom_vvl_sf( zsshvn_e, 'V', zfse3vn_e ) 
    535531 
    536532            ! Ocean depth at U- and V-points 

  • branches/dev_004_VVL/NEMO/OPA_SRC/TRA/trazdf.F90

    r1239 r1380  
    7878         ztrds(:,:,:) = sa(:,:,:) 
    7979      ENDIF 
    80  
    81       IF( lk_vvl )   CALL dom_vvl_ssh( kt )      ! compute ssha field 
    8280 
    8381      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend 

  • branches/dev_004_VVL/NEMO/OPA_SRC/istate.F90

    r1368 r1380  
    156156         CALL bn2( tb, sb, rn2 )    ! before Brunt Vaissala frequency 
    157157 
     158         IF( .NOT. ln_rstart ) CALL wzv( nit000 ) 
     159 
    158160      ENDIF 
     161 
     162      !                                       ! Vertical velocity 
     163      !                                       ! ----------------- 
     164 
     165      IF( .NOT. lk_vvl )    CALL wzv( nit000 )                         ! from horizontal divergence 
    159166 
    160167   END SUBROUTINE istate_init 

  • branches/dev_004_VVL/NEMO/OPA_SRC/oce.F90

    r1152 r1380  
    5959      gtv, gsv, grv          !: v-points at bottom ocean level  
    6060 
    61    !! free surface 
    62    !! ------------ 
    63    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !: 
    64       sshb , sshn ,        &  !: before, now sea surface height (meters) 
    65       sshu , sshv ,        &  !: sea surface height at u- and v- point 
    66       sshbb, ssha             !: before before sea surface height at t-point 
     61   !! free surface                       !  before  !  now     !  after   ! 
     62   !! ------------                       !  fields  !  fields  !  trends  ! 
     63   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshb   ,  sshn    ,  ssha    !: sea surface height at t-point [m] 
     64   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshu_b ,  sshu_n  ,  sshu_a  !: sea surface height at u-point [m] 
     65   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshv_b ,  sshv_n  ,  sshv_a  !: sea surface height at u-point [m] 
     66   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshf_b ,  sshf_n  ,  sshf_a  !: sea surface height at f-point [m] 
     67 
     68   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshbb         !: before before & after sea surface height at t-point 
     69   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshn_f        !: filtered now sea surface height at t-point 
    6770 
    6871#if defined key_dynspg_rl   ||   defined key_esopa 

  • branches/dev_004_VVL/NEMO/OPA_SRC/step.F90

    r1368 r1380  
    209209      IF( n_cla == 1 ) CALL div_cla( kstp )                 ! Cross Land Advection (Update Hor. divergence) 
    210210                       CALL ssh_wzv( kstp       )           ! after ssh & vertical velocity 
    211       IF( lk_vvl )     CALL dom_vvl                         ! vertical mesh at next time step 
    212  
    213211 
    214212      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
Note: See TracChangeset for help on using the changeset viewer.