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

Changeset 1389


Ignore:
Timestamp:
2009-04-07T17:54:50+02:00 (15 years ago)
Author:
rblod
Message:

Cleaning in VVL (initialisation), see ticket #402

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

Legend:

Unmodified
Added
Removed
  • branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domvvl.F90

    r1385 r1389  
    1212   !!---------------------------------------------------------------------- 
    1313   !!---------------------------------------------------------------------- 
    14    !!   dom_vvl         : defined scale factors & depths at each time step 
     14   !!   dom_vvl         : empty routine 
    1515   !!   dom_vvl_ini     : defined coefficients to distribute ssh on each layers 
    1616   !!---------------------------------------------------------------------- 
     
    192192 !     ENDIF 
    193193 
    194       ! Scale factors at T levels 
    195       DO jk = 1, jpkm1 
    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 
    214  
    215       ! Local depth or Inverse of the local depth of the water column at u- and v-points 
    216       ! ------------------------------ 
    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) 
     194       WRITE(*,*) 'dom_vvl  : empty routine, you should not be here' 
    223195 
    224196   END SUBROUTINE dom_vvl 
  • branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r1381 r1389  
    408408           CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) ) 
    409409           CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) ) 
    410            CALL iom_get( numror, jpdom_autoglo, 'sshb', sshb(:,:)         ) 
    411            CALL iom_get( numror, jpdom_autoglo, 'sshn', sshn(:,:)         ) 
    412410           IF( neuler == 0 ) THEN 
    413               sshb(:,:) = sshn(:,:) 
    414411              gcxb(:,:) = gcx (:,:) 
    415412           ENDIF 
     
    417414           gcx (:,:) = 0.e0 
    418415           gcxb(:,:) = 0.e0 
    419            IF( nn_rstssh == 1 ) THEN   
    420               sshb(:,:) = 0.e0 
    421               sshn(:,:) = 0.e0 
    422            ENDIF 
    423416        ENDIF 
    424417     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
     
    427420        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) ) 
    428421        CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) ) 
    429         CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:)         ) 
    430         CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:)         ) 
    431422     ENDIF 
    432423     ! 
  • branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r1380 r1389  
    607607      ! 
    608608      IF( TRIM(cdrw) == 'READ' ) THEN 
    609          IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 
    610             CALL iom_get( numror, jpdom_autoglo, 'sshb'  , sshb(:,:)   ) 
    611             CALL iom_get( numror, jpdom_autoglo, 'sshn'  , sshn(:,:)   ) 
    612             IF( neuler == 0 ) sshb(:,:) = sshn(:,:) 
    613          ELSE 
    614             IF( nn_rstssh == 1 ) THEN   
    615                sshb(:,:) = 0.e0 
    616                sshn(:,:) = 0.e0 
    617             ENDIF 
    618          ENDIF 
    619609         IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN 
    620610            CALL iom_get( numror, jpdom_autoglo, 'sshb_b', sshb_b(:,:) )   ! free surface issued 
     
    644634         ENDIF 
    645635      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
    646          CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb  (:,:) ) 
    647          CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn  (:,:) ) 
    648636         CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) )   ! free surface issued 
    649637         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) )   ! from barotropic loop 
  • branches/dev_004_VVL/NEMO/OPA_SRC/DYN/wzvmod.F90

    r1386 r1389  
    99   !!              "   !  07-07  (D. Storkey) Zero zhdiv at open boundary (BDY)  
    1010   !!---------------------------------------------------------------------- 
    11    !!   wzv        : Compute the vertical velocity 
     11   !!   wzv            : empty routine  
     12   !!   ssh_wzv        : after ssh & now vertical velocity 
     13   !!   ssh_nxt        : filter ans swap the ssh arrays 
     14   !!   ssh_rst        : read/write ssh restart fields in the ocean restart file 
    1215   !!---------------------------------------------------------------------- 
    1316   !! * Modules used 
     
    1619   USE sbc_oce         ! surface boundary condition: ocean 
    1720   USE domvvl          ! Variable volume 
     21   USE iom             ! I/O library 
     22   USE restart         ! only for lrst_oce 
    1823   USE in_out_manager  ! I/O manager 
    1924   USE prtctl          ! Print control 
     
    5964      !!---------------------------------------------------------------------- 
    6065      INTEGER, INTENT(in) :: kt 
    61       !! 
    62       INTEGER  ::   jk     ! dummy loop indices 
    63       REAL(wp) ::   z2dt   ! temporary scalar 
    64       !!---------------------------------------------------------------------- 
    65  
    66       IF( kt == nit000 ) THEN 
    67          IF(lwp) WRITE(numout,*) 
    68          IF(lwp) WRITE(numout,*) 'wzv     : vertical velocity from continuity eq.' 
    69          IF(lwp) WRITE(numout,*) '~~~~~~~ '  
    70  
    71          ! bottom boundary condition: w=0 (set once for all) 
    72          wn(:,:,jpk) = 0.e0 
    73       ENDIF 
    74  
    75       IF( lk_vvl ) THEN                ! Variable volume 
    76          ! 
    77          z2dt = 2. * rdt                                       ! time step: leap-frog 
    78          IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest 
    79  
    80  
    81          ! Vertical velocity computed from bottom 
    82          ! -------------------------------------- 
    83          DO jk = jpkm1, 1, -1 
    84             wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) & 
    85               &        - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) / z2dt 
    86          END DO 
    87  
    88       ELSE                             ! Fixed volume 
    89  
    90          ! Vertical velocity computed from bottom 
    91          ! -------------------------------------- 
    92          DO jk = jpkm1, 1, -1 
    93             wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) 
    94          END DO 
    95  
    96       ENDIF 
    97  
    98       IF(ln_ctl)   CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 -   : ', mask1=wn) 
    99       ! 
     66       
     67      ! Empty routine 
     68 
     69      WRITE(*,*) 'wzv : you should not be here : error ?'   
     70 
    10071   END SUBROUTINE wzv 
    10172 
     
    10374   SUBROUTINE ssh_wzv( kt )  
    10475      !!---------------------------------------------------------------------- 
    105       !!                ***  ROUTINE dom_vvl_ssh  *** 
     76      !!                ***  ROUTINE ssh_wzv  *** 
    10677      !!                    
    10778      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity 
    108       !!              Cand update the now vertical coordinate (lk_vvl=T). 
    109       !! 
    110       !! ** Method  : -   
    111  
    112       !!              - Using the incompressibility hypothesis, the vertical 
    113       !!      velocity is computed by integrating the horizontal divergence  
     79      !!              and update the now vertical coordinate (lk_vvl=T). 
     80      !! 
     81      !! ** Method  : - 
     82      !! 
     83      !!              - Using the incompressibility hypothesis, the vertical  
     84      !!      velocity is computed by integrating the horizontal divergence   
    11485      !!      from the bottom to the surface minus the scale factor evolution. 
    11586      !!        The boundary conditions are w=0 at the bottom (no flux) and. 
    11687      !! 
    117       !! ** action  :   wn array : the now vertical velocity 
     88      !! ** action  :   ssha    : after sea surface height 
     89      !!                wn      : now vertical velocity 
     90      !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height 
     91      !!                          at u-, v-, f-point s 
     92      !!                .._1    : now vertical coordinate arrays 
     93      !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 
    11894      !!---------------------------------------------------------------------- 
    11995      INTEGER, INTENT(in) ::   kt   ! time step 
    12096      !! 
    12197      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
     98      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars 
    12299      REAL(wp) ::   z2dt, zraur     ! temporary scalars 
    123100      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace 
     
    129106         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    130107         ! 
    131          wn(:,:,jpk) = 0.e0      ! bottom boundary condition: w=0 (set once for all) 
     108         CALL ssh_rst( nit000, 'READ' )       ! read or initialize sshb and sshn 
     109         ! 
     110         wn(:,:,jpk) = 0.e0                   ! bottom boundary condition: w=0 (set once for all) 
     111         ! 
     112         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only) 
     113            DO jj = 1, jpjm1 
     114               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible 
     115                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
     116                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
     117                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 
     118                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
     119                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
     120                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
     121                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
     122                  sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 & 
     123                     &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) ) 
     124                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
     125                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
     126                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
     127                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
     128                  sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 & 
     129                     &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) 
     130               END DO 
     131            END DO 
     132            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
     133            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
     134            CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. ) 
     135         ENDIF 
     136         ! 
    132137      ENDIF 
    133138 
     
    173178      !                                           !------------------------------! 
    174179      !                                                ! integrate from the bottom the hor. divergence 
    175          DO jk = jpkm1, 1, -1 
    176             wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   & 
    177                  &                    - (  fse3t_a(:,:,jk)                   & 
    178                  &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 
    179          END DO 
     180      DO jk = jpkm1, 1, -1 
     181         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   & 
     182              &                    - (  fse3t_a(:,:,jk)                   & 
     183              &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 
     184      END DO 
    180185 
    181186      !                                           !------------------------------! 
     
    201206         hv(:,:) = hv_0(:,:) + sshv_n(:,:) 
    202207         !                                             ! masked inverse of the ocean depth (at u- and v-points) 
     208         hur(:,:) = 1. / MAX( hu(:,:), fse3u_0(:,:,1) ) * umask(:,:,1) 
     209         hvr(:,:) = 1. / MAX( hv(:,:), fse3v_0(:,:,1) ) * vmask(:,:,1) 
     210!!gm to be corrected (the above case does not work properly with 1 ocean level only) 
     211!        hur(:,:) = 1. / MAX( hu(:,:), 1.e-15 ) * umask(:,:,1) 
     212!        hvr(:,:) = 1. / MAX( hv(:,:), 1.e-15 ) * vmask(:,:,1) 
     213!!gm 
     214!!add end 
     215 
    203216      ENDIF 
    204217      ! 
     
    273286      ENDIF 
    274287 
    275  
     288      !                      ! write filtered free surface arrays in restart file 
     289      IF( lrst_oce )   CALL ssh_rst( kt, 'WRITE' ) 
     290      ! 
    276291      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 ) 
    277292      ! 
    278293   END SUBROUTINE ssh_nxt 
     294 
     295 
     296   SUBROUTINE ssh_rst( kt, cdrw ) 
     297      !!--------------------------------------------------------------------- 
     298      !!                   ***  ROUTINE ssh_rst  *** 
     299      !! 
     300      !! ** Purpose :   Read or write Sea Surface Height arrays in restart file 
     301      !! 
     302      !! ** action  : - cdrw = READ  :  sshb, sshn  read    in ocean restart file 
     303      !!              - cdrw = WRITE :  sshb, sshn  written in ocean restart file 
     304      !!---------------------------------------------------------------------- 
     305      INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
     306      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
     307      !!---------------------------------------------------------------------- 
     308      ! 
     309      IF( TRIM(cdrw) == 'READ' ) THEN 
     310         IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 
     311            CALL iom_get( numror, jpdom_autoglo, 'sshb'  , sshb(:,:)   ) 
     312            CALL iom_get( numror, jpdom_autoglo, 'sshn'  , sshn(:,:)   ) 
     313            IF( neuler == 0 )   sshb(:,:) = sshn(:,:) 
     314         ELSE 
     315            IF( nn_rstssh == 1 ) THEN 
     316               sshb(:,:) = 0.e0 
     317               sshn(:,:) = 0.e0 
     318            ENDIF 
     319         ENDIF 
     320      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
     321         CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb  (:,:) ) 
     322         CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn  (:,:) ) 
     323      ENDIF 
     324      ! 
     325   END SUBROUTINE ssh_rst 
    279326 
    280327   !!====================================================================== 
  • branches/dev_004_VVL/NEMO/OPA_SRC/istate.F90

    r1380 r1389  
    132132          
    133133      ENDIF 
    134  
    135       IF( lk_vvl .OR. lk_agrif ) THEN 
     134      ! 
     135      IF( lk_agrif ) THEN 
    136136         ! read free surface arrays in restart file 
    137137         IF( ln_rstart ) THEN 
    138138            IF( lk_dynspg_flt )   CALL flt_rst( nit000, 'READ' )      ! read or initialize the following fields 
    139             !                                                         ! gcx, gcxb, sshb, sshn 
    140             IF( lk_dynspg_ts  )   CALL ts_rst ( nit000, 'READ' )      ! read or initialize the following fields 
    141             !                                                         ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 
    142             IF( lk_dynspg_exp )   CALL exp_rst( nit000, 'READ' )      ! read or initialize the following fields 
    143             !                                                         ! sshb, sshn 
    144          ENDIF 
     139            !                                                         ! gcx, gcxb for agrif_opa_init 
     140         ENDIF                                                        ! explicit case not coded yet with AGRIF 
    145141      ENDIF 
    146  
    147       IF( lk_vvl ) THEN 
    148  
    149          ! 
    150          IF( .NOT. lk_dynspg_flt ) sshbb(:,:) = sshb(:,:) 
    151          ! 
    152          CALL dom_vvl               ! ssh init and vertical grid update 
    153  
    154          CALL eos_init              ! usefull to get the equation state type neos parameter 
    155  
    156          CALL bn2( tb, sb, rn2 )    ! before Brunt Vaissala frequency 
    157  
    158          IF( .NOT. ln_rstart ) CALL wzv( nit000 ) 
    159  
    160       ENDIF 
    161  
    162       !                                       ! Vertical velocity 
    163       !                                       ! ----------------- 
    164  
    165       IF( .NOT. lk_vvl )    CALL wzv( nit000 )                         ! from horizontal divergence 
    166  
     142      ! 
    167143   END SUBROUTINE istate_init 
    168144 
Note: See TracChangeset for help on using the changeset viewer.