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

Changeset 2068


Ignore:
Timestamp:
2010-09-06T17:56:51+02:00 (14 years ago)
Author:
mlelod
Message:

ticket: #663 ensuring restartability and conservation

Location:
branches/DEV_r1837_MLF/NEMO/OPA_SRC
Files:
14 edited

Legend:

Unmodified
Added
Removed
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/dom_oce.F90

    r1792 r2068  
    145145   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3vw_1             !: analytical vertical scale factors at  VW-- 
    146146   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3w_1  , e3uw_1    !:                                       W--UW  points (m) 
     147   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3t_b              !: before         -      -      -    -   T      points (m) 
     148   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3u_b  , e3v_b     !:   -            -      -      -    -   U--V   points (m) 
     149   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3t_d              !: "diffused"     -      -      -    -   T      points (m) 
    147150#else 
    148151   LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domvvl.F90

    r2005 r2068  
    130130            zcoeff = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj) ) 
    131131            ! before fields 
    132             zv_t_ij       = zs_t(ji  ,jj  ) * sshb(ji  ,jj  ,jk) 
    133             zv_t_ip1j     = zs_t(ji+1,jj  ) * sshb(ji+1,jj  ,jk) 
    134             zv_t_ijp1     = zs_t(ji  ,jj+1) * sshb(ji  ,jj+1,jk) 
     132            zv_t_ij       = zs_t(ji  ,jj  ) * sshb(ji  ,jj  ) 
     133            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshb(ji+1,jj  ) 
     134            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshb(ji  ,jj+1) 
    135135            sshu_b(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j ) 
    136136            sshv_b(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 ) 
    137137            ! now fields 
    138             zv_t_ij       = zs_t(ji  ,jj  ) * sshn(ji  ,jj  ,jk) 
    139             zv_t_ip1j     = zs_t(ji+1,jj  ) * sshn(ji+1,jj  ,jk) 
    140             zv_t_ijp1     = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1,jk) 
    141             zv_t_ip1jp1   = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1,jk) 
     138            zv_t_ij       = zs_t(ji  ,jj  ) * sshn(ji  ,jj  ) 
     139            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshn(ji+1,jj  ) 
     140            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1) 
     141            zv_t_ip1jp1   = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1) 
    142142            sshu_n(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j ) 
    143143            sshv_n(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 ) 
     
    162162         END DO 
    163163      END DO 
    164       CALL lbc_lnk( fse3u_b, 'U', 1. )               ! lateral boundary conditions 
    165       CALL lbc_lnk( fse3v_b, 'U', 1. ) 
     164      CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions 
     165      CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 
    166166      ! Add initial scale factor to scale factor anomaly 
    167167      fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/DOM/domzgr_substitute.h90

    r2005 r2068  
    6363#   define  fse3vw_n(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_n(i,j)*muv(i,j,k))) 
    6464 
    65 #   define  fse3t_m(i,j,k)   e3t_m(i,j,k) 
     65#   define  fse3t_m(i,j,k)   (fse3t_0(i,j,k)*(1.+ssh_m(i,j)*mut(i,j,k))) 
     66#   define  fse3t_d(i,j,k)   e3t_d(i,j,k) 
    6667 
    6768#   define  fse3t_a(i,j,k)   (fse3t_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 
     
    100101 
    101102#   define  fse3t_m(i,j,k)   fse3t_0(i,j,k) 
     103#   define  fse3t_d(i,j,k)   fse3t_0(i,j,k) 
    102104 
    103105#   define  fse3t_a(i,j,k)   fse3t_0(i,j,k) 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90

    r2005 r2068  
    2222   USE oce             ! ocean dynamics and tracers 
    2323   USE dom_oce         ! ocean space and time domain 
     24   USE sbc_oce         ! Surface boundary condition: ocean fields 
     25   USE phycst          ! physical constants 
    2426   USE dynspg_oce      ! type of surface pressure gradient 
    2527   USE dynadv          ! dynamics: vector invariant versus flux form 
     
    218220         END DO 
    219221      ELSE                                             !* Leap-Frog : Asselin filter and swap 
    220          IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN          ! applied on velocity 
     222         !                                ! =============! 
     223         IF( .NOT. lk_vvl ) THEN          ! Fixed volume ! 
     224            !                             ! =============! 
    221225            DO jk = 1, jpkm1                               
    222226               DO jj = 1, jpj 
     
    232236               END DO 
    233237            END DO 
    234          ELSE                                                ! applied on thickness weighted velocity 
    235             ! Before scale factors at (t/u/v)-points (actually "now filtered" and futur "before") 
    236             ! ====================================== 
    237             ! Scale factor at t-points 
    238             ! ------------------------ 
    239             fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * fse3t_m(:,:,:) 
     238            !                             ! ================! 
     239         ELSE                             ! Variable volume ! 
     240            !                             ! ================! 
     241            ! Before scale factor at t-points 
     242            ! ------------------------------- 
     243            DO jk = 1, jpkm1 
     244               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) + atfp * fse3t_d(:,:,jk) 
     245            ENDDO 
    240246            ! Add volume filter correction only at the first level of t-point scale factors 
    241247            zec = atfp * rdt / rau0 
    242248            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
    243             ! Scale factor at (u/v)-points 
    244             ! ------------------------ 
    245249            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 
    246250            zs_t  (:,:) =       e1t(:,:) * e2t(:,:) 
    247251            zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:) 
    248252            zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:) 
    249             ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
    250             DO jk = 1, jpkm1 
    251                DO jj = 1, jpjm1 
    252                   DO ji = 1, jpim1 
    253                      zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
    254                      zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
    255                      zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
    256                      ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
    257                      ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
    258                   END DO 
    259                END DO 
    260             END DO 
    261             CALL lbc_lnk( ze3u_f, 'U', 1. )               ! lateral boundary conditions 
    262             CALL lbc_lnk( ze3v_f, 'U', 1. ) 
    263             ! Add initial scale factor to scale factor anomaly 
    264             ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:) 
    265             ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:) 
    266              
    267             DO jk = 1, jpkm1 
    268                DO jj = 1, jpj 
    269                   DO ji = 1, jpim 
    270                      zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk) 
    271                      zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk) 
    272                      zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk) 
    273                      zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk) 
    274                      zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk) 
    275                      zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk) 
    276                      ! 
    277                      zuf  = ( zue3n + atfp * ( zue3b  - 2.e0 * zue3n  + zue3a  ) ) / ze3u_f(ji,jj,jk) 
    278                      zvf  = ( zve3n + atfp * ( zve3b  - 2.e0 * zve3n  + zve3a  ) ) / ze3v_f(ji,jj,jk) 
    279                      ! 
    280                      ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
    281                      vb(ji,jj,jk) = zvf 
    282                      fse3u_b(ji,jj,jk) = ze3u_f(ji,jj,jk)    ! e3u_b <-- filtered scale factor 
    283                      fse3v_b(ji,jj,jk) = ze3v_f(ji,jj,jk) 
    284                      un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
    285                      vn(ji,jj,jk) = va(ji,jj,jk) 
    286                   END DO 
    287                END DO 
    288             END DO 
    289             CALL lbc_lnk( ub, 'U', -1. )         ! local domain boundaries 
    290             CALL lbc_lnk( vb, 'V', -1. )  
     253            ! 
     254            IF( ln_dynadv_vec ) THEN 
     255               ! Before scale factor at (u/v)-points 
     256               ! ----------------------------------- 
     257               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
     258               DO jk = 1, jpkm1 
     259                  DO jj = 1, jpjm1 
     260                     DO ji = 1, jpim1 
     261                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
     262                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
     263                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
     264                        fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
     265                        fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
     266                     END DO 
     267                  END DO 
     268               END DO 
     269               ! lateral boundary conditions 
     270               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. ) 
     271               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 
     272               ! Add initial scale factor to scale factor anomaly 
     273               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 
     274               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 
     275               ! Leap-Frog - Asselin filter and swap: applied on velocity 
     276               ! ----------------------------------- 
     277               DO jk = 1, jpkm1 
     278                  DO jj = 1, jpj 
     279                     DO ji = 1, jpi 
     280                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) ) 
     281                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) ) 
     282                        ! 
     283                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
     284                        vb(ji,jj,jk) = zvf 
     285                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
     286                        vn(ji,jj,jk) = va(ji,jj,jk) 
     287                     END DO 
     288                  END DO 
     289               END DO 
     290               ! 
     291            ELSE 
     292               ! Temporary filered scale factor at (u/v)-points (will become before scale factor) 
     293               !----------------------------------------------- 
     294               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
     295               DO jk = 1, jpkm1 
     296                  DO jj = 1, jpjm1 
     297                     DO ji = 1, jpim1 
     298                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
     299                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
     300                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
     301                        ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
     302                        ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
     303                     END DO 
     304                  END DO 
     305               END DO 
     306               ! lateral boundary conditions 
     307               CALL lbc_lnk( ze3u_f, 'U', 1. ) 
     308               CALL lbc_lnk( ze3v_f, 'V', 1. ) 
     309               ! Add initial scale factor to scale factor anomaly 
     310               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:) 
     311               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:) 
     312               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 
     313               ! -----------------------------------             =========================== 
     314               DO jk = 1, jpkm1 
     315                  DO jj = 1, jpj 
     316                     DO ji = 1, jpim1 
     317                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk) 
     318                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk) 
     319                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk) 
     320                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk) 
     321                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk) 
     322                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk) 
     323                        ! 
     324                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
     325                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
     326                        ! 
     327                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
     328                        vb(ji,jj,jk) = zvf 
     329                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
     330                        vn(ji,jj,jk) = va(ji,jj,jk) 
     331                     END DO 
     332                  END DO 
     333               END DO 
     334               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor 
     335               fse3v_b(:,:,:) = ze3v_f(:,:,:) 
     336               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions 
     337               CALL lbc_lnk( vb, 'V', -1. ) 
     338            ENDIF 
     339            ! 
    291340         ENDIF 
     341         ! 
    292342      ENDIF 
    293343 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90

    r2005 r2068  
    104104            DO jj = 1, jpjm1 
    105105               DO ji = 1, jpim1      ! NO Vector Opt. 
    106                   sshf_b(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
    107                      &                 / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
    108                      &                 * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_b(ji,jj  )     & 
    109                      &                   + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_b(ji,jj+1) ) 
    110                END DO 
    111             END DO 
    112             DO jj = 1, jpjm1 
    113                DO ji = 1, jpim1      ! NO Vector Opt. 
    114106                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
    115107                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
     
    118110               END DO 
    119111            END DO 
    120             CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. ) 
     112            CALL lbc_lnk( sshf_n, 'F', 1. ) 
    121113         ENDIF 
    122114         ! 
     
    157149      !                                           !   After Sea Surface Height   ! 
    158150      !                                           !------------------------------! 
    159       z1_rau0 = 0.5 / rau0 
    160       ! 
    161151      zhdiv(:,:) = 0.e0 
    162152      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
     
    166156      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used 
    167157      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp 
    168       ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     158      IF( neuler == 0 .AND. kt == nit000 ) THEN 
     159         z1_rau0 = 1.e0 / rau0 
     160         ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 *   emp(:,:)                + zhdiv(:,:) )  ) * tmask(:,:,1) 
     161      ELSE 
     162         z1_rau0 = 0.5 / rau0 
     163         ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     164      ENDIF 
    169165 
    170166#if defined key_obc 
     
    194190      !                                           !----------------------------------------! 
    195191      ! Needed for Robert-Asselin time filter and for Brown & Campana semi implicit hydrostatic presure gradient 
    196       fse3t_m(:,:,:) =          fse3t_b(:,:,:)   & 
    197          &             - 2.e0 * fse3t_n(:,:,:)   & 
    198          &             +        fse3t_a(:,:,:) 
     192      DO jk = 1, jpk 
     193         fse3t_d(:,:,jk) =          fse3t_b(:,:,jk)   & 
     194            &              - 2.e0 * fse3t_n(:,:,jk)   & 
     195            &              +        fse3t_a(:,:,jk) 
     196      ENDDO 
    199197      !                                           !------------------------------! 
    200198      !                                           !     Now Vertical Velocity    ! 
     
    251249      !! 
    252250      INTEGER  ::   ji, jj   ! dummy loop indices 
    253       REAL(wp) ::   zec      ! temporary scalare 
     251      REAL(wp) ::   zec      ! temporary scalar 
    254252      !!---------------------------------------------------------------------- 
    255253 
     
    281279            CALL lbc_lnk( sshf_n, 'F', 1. ) 
    282280         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap 
     281            zec = atfp * rdt / rau0 
    283282            DO jj = 1, jpj 
    284283               DO ji = 1, jpi                                ! before <-- now filtered 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90

    r2005 r2068  
    2626   USE zdfmxl          ! mixed layer depth 
    2727   USE trdmld_oce      ! ocean active mixed layer tracers trends variables 
     28   USE domvvl          ! variable volume 
    2829 
    2930   IMPLICIT NONE 
     
    3839 
    3940   !! * Substitutions 
     41#  include "domzgr_substitute.h90" 
    4042#  include "vectopt_loop_substitute.h90" 
    4143   !!---------------------------------------------------------------------- 
     
    123125      CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb    ) 
    124126      IF( lk_vvl ) THEN 
    125          CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b ) 
     127         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) ) 
    126128      ENDIF 
    127129      ! 
     
    157159      REAL(wp) ::   zrdt, zrdttra1 
    158160      INTEGER  ::   jlibalt = jprstlib 
     161      INTEGER  ::   jk                    ! dummy loop indices 
    159162      LOGICAL  ::   llok 
    160163      !!---------------------------------------------------------------------- 
     
    195198      CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb  ) 
    196199      IF( lk_vvl ) THEN 
    197          CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b ) 
     200         CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
    198201      ENDIF 
    199202      ! 
     
    224227         sshb (:,:)   = sshn (:,:) 
    225228         IF( lk_vvl ) THEN 
    226             fse3t_b(:,:,:) = fse3t_n(:,:,:) 
     229            DO jk = 1, jpk 
     230               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 
     231            ENDDO 
    227232         ENDIF 
    228233      ENDIF 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r1975 r2068  
    5252   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps   , emps_b   !: freshwater budget: concentration/dillution   [Kg/m2/s] 
    5353   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot           !: total E-P-R over ocean and ice               [Kg/m2/s] 
    54    ! - ML - beginning 
     54   ! - ML - begin 
    5555   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_trd_hc_n      !: sbc heat content trend now                   [K/m/s] 
    5656   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_trd_hc_b      !:  "   "      "      "   before                   " 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90

    r1975 r2068  
    194194         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields 
    195195         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields 
    196          qns_b (:,:) = qns (:,:)                         !  are set the end of the routine) 
     196         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine) 
    197197         qsr_b (:,:) = qsr (:,:) 
    198198         emp_b (:,:) = emp (:,:) 
    199199         emps_b(:,:) = emps(:,:) 
    200200         ! - ML - 
    201          sbc_trd_hc_b(:,:) = sbc_trd_hc_n(:,:) 
    202          qsr_trd_hc_b(:,:,:) = qsr_trd_hc_n(:,:,:) 
    203          IF ( .NOT. lk_vvl )  sbc_trd_sc_b(:,:) = sbc_trd_sc_n(:,:) 
     201         ! sbc_trd_hc_b(:,:) = sbc_trd_hc_n(:,:) 
     202         ! qsr_trd_hc_b(:,:,:) = qsr_trd_hc_n(:,:,:) 
     203         ! IF ( .NOT. lk_vvl )  sbc_trd_sc_b(:,:) = sbc_trd_sc_n(:,:) 
    204204          
    205205      ENDIF 
     
    265265            CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  )   ! before     solar heat flux (T-point) 
    266266            CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b  )   ! before     freshwater flux (T-point) 
    267             CALL iom_get( numror, jpdom_autoglo, 'emps_b', emp_b )   ! before C/D freshwater flux (T-point) 
     267            CALL iom_get( numror, jpdom_autoglo, 'emps_b', emps_b )   ! before C/D freshwater flux (T-point) 
    268268            ! - ML - 
    269             CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_hc_b', sbc_trd_hc_b )   ! before heat content sbc trend 
    270             CALL iom_get( numror, jpdom_autoglo, 'qsr_trd_hc_b', qsr_trd_hc_b )   ! before heat content trend due to Qsr flux 
    271             IF ( .NOT. lk_vvl ) THEN 
    272                CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_sc_b', sbc_trd_sc_b )   ! before salt content sbc trend 
    273             ENDIF 
     269            ! CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_hc_b', sbc_trd_hc_b )   ! before heat content sbc trend 
     270            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_trd_hc_b', qsr_trd_hc_b )   ! before heat content trend due to Qsr flux 
     271            ! IF ( .NOT. lk_vvl ) THEN 
     272            !    CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_sc_b', sbc_trd_sc_b )   ! before salt content sbc trend 
     273            ! ENDIF 
    274274            ! 
    275275         ELSE                                                   !* no restart: set from nit000 values 
     
    291291            &                    'at it= ', kt,' date= ', ndastp 
    292292         IF(lwp) WRITE(numout,*) '~~~~' 
    293          CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau_b )    !  
    294          CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau_b ) 
    295          CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns_b  ) 
    296          CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr_b  ) 
    297          CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp_b  ) 
    298          CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps_b ) 
     293         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )    !  
     294         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau ) 
     295         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  ) 
     296         CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  ) 
     297         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  ) 
     298         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps ) 
    299299         ! - ML - 
    300          CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_hc_b', sbc_trd_hc_b ) 
    301          CALL iom_rstput( kt, nitrst, numrow, 'qsr_trd_hc_b', qsr_trd_hc_b ) 
    302          IF ( .NOT. lk_vvl ) THEN 
    303             CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_sc_b', sbc_trd_sc_b ) 
    304          ENDIF 
     300         ! CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_hc_b', sbc_trd_hc_n ) 
     301         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_trd_hc_b', qsr_trd_hc_n ) 
     302         ! IF ( .NOT. lk_vvl ) THEN 
     303         !    CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_sc_b', sbc_trd_sc_n ) 
     304         ! ENDIF 
    305305         ! 
    306306      ENDIF 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trabbc.F90

    r1601 r2068  
    3535   REAL(wp) ::   rn_geoflx_cst = 86.4e-3      ! Constant value of geothermal heat flux 
    3636 
    37    INTEGER , DIMENSION(jpi,jpj) ::   nbotlevt   ! ocean bottom level index at T-pt 
    38    REAL(wp), DIMENSION(jpi,jpj) ::   qgh_trd0   ! geothermal heating trend 
     37   INTEGER , PUBLIC, DIMENSION(jpi,jpj) ::   nbotlevt   ! ocean bottom level index at T-pt 
     38   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qgh_trd0   ! geothermal heating trend 
    3939  
    4040   !! * Substitutions 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90

    r2005 r2068  
    178178      !! 
    179179      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
    180       REAL(wp) ::   zt_m, zs_m      ! temporary scalars 
     180      REAL(wp) ::   zt_d, zs_d      ! temporary scalars 
    181181      REAL(wp) ::   ztn, zsn        !    -         - 
    182182      !!---------------------------------------------------------------------- 
     
    203203                  !                                         ! time laplacian on tracers 
    204204                  !                                         ! used for both Asselin and Brown & Campana filters 
    205                   zt_m = ta(ji,jj,jk) - 2. * tn(ji,jj,jk) + tb(ji,jj,jk) 
    206                   zs_m = sa(ji,jj,jk) - 2. * sn(ji,jj,jk) + sb(ji,jj,jk) 
     205                  zt_d = ta(ji,jj,jk) - 2. * tn(ji,jj,jk) + tb(ji,jj,jk) 
     206                  zs_d = sa(ji,jj,jk) - 2. * sn(ji,jj,jk) + sb(ji,jj,jk) 
    207207                  ! 
    208208                  !                                         ! swap of arrays 
    209                   tb(ji,jj,jk) = tn(ji,jj,jk) + atfp * zt_m               ! tb <-- tn filtered 
    210                   sb(ji,jj,jk) = sn(ji,jj,jk) + atfp * zs_m               ! sb <-- sn filtered 
     209                  tb(ji,jj,jk) = tn(ji,jj,jk) + atfp * zt_d               ! tb <-- tn filtered 
     210                  sb(ji,jj,jk) = sn(ji,jj,jk) + atfp * zs_d               ! sb <-- sn filtered 
    211211                  tn(ji,jj,jk) = ta(ji,jj,jk)                             ! tn <-- ta 
    212212                  sn(ji,jj,jk) = sa(ji,jj,jk)                             ! sn <-- sa 
    213213                  !                                         ! semi imlicit hpg computation (Brown & Campana) 
    214214                  IF( ln_dynhpg_imp ) THEN 
    215                      ta(ji,jj,jk) = ztn + rbcp * zt_m                     ! ta <-- Brown & Campana average 
    216                      sa(ji,jj,jk) = zsn + rbcp * zs_m                     ! sa <-- Brown & Campana average 
     215                     ta(ji,jj,jk) = ztn + rbcp * zt_d                     ! ta <-- Brown & Campana average 
     216                     sa(ji,jj,jk) = zsn + rbcp * zs_d                     ! sa <-- Brown & Campana average 
    217217                  ENDIF 
    218218               END DO 
     
    254254      REAL     ::   ztc_a, ztc_n, ztc_b            !    -         - 
    255255      REAL     ::   zsc_a, zsc_n, zsc_b            !    -         - 
    256       REAL     ::   ztc_f, zsc_f, ztc_m, zsc_m     !    -         - 
    257       REAL     ::   ze3t_f, ze3t_m                 !    -         - 
     256      REAL     ::   ztc_f, zsc_f, ztc_d, zsc_d     !    -         - 
     257      REAL     ::   ze3t_f, ze3t_d                 !    -         - 
    258258      REAL     ::   zfact1, zfact2                 !    -         - 
    259259      !!---------------------------------------------------------------------- 
     
    274274      ELSE                                             ! apply filter on thickness weighted tracer and swap 
    275275         DO jk = 1, jpkm1 
    276             zfact1 = atfp * r2dt_t(jk) 
     276            zfact1 = atfp * rdttra(jk) 
    277277            zfact2 = zfact1 / rau0 
    278278            DO jj = 1, jpj 
     
    282282                  ze3t_n = fse3t_n(ji,jj,jk) 
    283283                  ze3t_a = fse3t_a(ji,jj,jk) 
    284                   ze3t_m = fse3t_m(ji,jj,jk) 
     284                  ze3t_d = fse3t_d(ji,jj,jk) 
    285285                  !                                         ! tracer content at Before, now and after 
    286286                  ztc_b  = tb(ji,jj,jk) * ze3t_b   ;   zsc_b = sb(ji,jj,jk) * ze3t_b 
     
    290290                  !                                         ! Time laplacian on tracer contents 
    291291                  !                                         ! used for both Asselin and Brown & Campana filters 
    292                   ztc_m  = ztc_a  - 2. * ztc_n + ztc_b 
    293                   zsc_m  = zsc_a  - 2. * zsc_n + zsc_b 
     292                  ztc_d  = ztc_a - 2. * ztc_n + ztc_b 
     293                  zsc_d  = zsc_a - 2. * zsc_n + zsc_b 
    294294                  !                                         ! Asselin Filter on thicknesses and tracer contents 
    295                   ze3t_f = ze3t_n + atfp * ze3t_m 
    296                   ztc_f  = ztc_n  + atfp * ztc_m 
    297                   zsc_f  = zsc_n  + atfp * zsc_m 
     295                  ze3t_f = ze3t_n + atfp * ze3t_d 
     296                  ztc_f  = ztc_n  + atfp * ztc_d 
     297                  zsc_f  = zsc_n  + atfp * zsc_d 
    298298                  !                                         ! Filter correction 
    299299                  IF( jk == 1 ) THEN 
    300                      ze3t_f = ze3t_f - zfact2 * ( emp_b       (ji,jj)    - emp         (ji,jj)    ) 
    301                      ztc_f  = ztc_f  - zfact1 * ( sbc_trd_hc_n(ji,jj)    - sbc_trd_hc_b(ji,jj)    ) 
     300                     ! WRITE(numout,*) 'filter correction: sbc_trd_hc_n' 
     301                     ze3t_f = ze3t_f - zfact2 * ( emp_b       (ji,jj) - emp         (ji,jj) ) 
     302                     ztc_f  = ztc_f  - zfact1 * ( sbc_trd_hc_n(ji,jj) - sbc_trd_hc_b(ji,jj) ) 
    302303                  ENDIF 
    303304                  IF( ln_traqsr .AND. ( jk .LE. nksr ) ) THEN 
     305                     ! WRITE(numout,*) 'jk =', jk 
     306                     ! WRITE(numout,*) 'filter correction: qsr_trd_hc_n' 
    304307                     ztc_f  = ztc_f  - zfact1 * ( qsr_trd_hc_n(ji,jj,jk) - qsr_trd_hc_b(ji,jj,jk) ) 
    305308                  ENDIF 
    306                   !                                         ! swap of arrays 
     309                                                          ! swap of arrays 
    307310                  ze3t_f = 1.e0 / ze3t_f 
    308311                  tb(ji,jj,jk) = ztc_f * ze3t_f                           ! tb <-- tn filtered 
     
    312315                  !                                         ! semi imlicit hpg computation (Brown & Campana) 
    313316                  IF( ln_dynhpg_imp ) THEN 
    314                      ze3t_m       = 1.e0   / ( ze3t_n + rbcp * ze3t_m ) 
    315                      ta(ji,jj,jk) = ze3t_m * ( ztc_n  + rbcp * ztc_m  )   ! ta <-- Brown & Campana average 
    316                      sa(ji,jj,jk) = ze3t_m * ( zsc_n  + rbcp * zsc_m  )   ! sa <-- Brown & Campana average 
     317                     ze3t_d       = 1.e0   / ( ze3t_n + rbcp * ze3t_d ) 
     318                     ta(ji,jj,jk) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average 
     319                     sa(ji,jj,jk) = ze3t_d * ( zsc_n  + rbcp * zsc_d  )   ! sa <-- Brown & Campana average 
    317320                  ENDIF 
    318321               END DO 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90

    r1975 r2068  
    112112         ztrdt(:,:,:) = ta(:,:,:)  
    113113         ztrds(:,:,:) = 0.e0 
     114      ENDIF 
     115 
     116      !                                           ! ---------------------------------------- ! 
     117      !                                           !          Swap of forcing field           ! 
     118      !                                           ! ---------------------------------------- ! 
     119      IF( kt /= nit000 ) qsr_trd_hc_b(:,:,:) = qsr_trd_hc_n(:,:,:) 
     120      !                                           ! ---------------------------------------- ! 
     121      IF( kt == nit000 ) THEN                     !   set the forcing field at nit000 - 1    ! 
     122         !                                        ! ---------------------------------------- ! 
     123         IF( ln_rstart .AND.    &                               !* Restart: read in restart file 
     124            & iom_varid( numror, 'qsr_trd_hc_b', ldstop = .FALSE. ) > 0 ) THEN  
     125            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file' 
     126            CALL iom_get( numror, jpdom_autoglo, 'qsr_trd_hc_b', qsr_trd_hc_b )   ! before heat content trend due to Qsr flux 
     127         ENDIF 
    114128      ENDIF 
    115129 
     
    225239            END DO 
    226240         END DO 
     241      ENDIF 
     242 
     243      !                                            ! ---------------------------------------- ! 
     244      IF( lrst_oce ) THEN                          !      Write in the ocean restart file     ! 
     245         !                                         ! ---------------------------------------- ! 
     246         IF(lwp) WRITE(numout,*) 
     247         IF(lwp) WRITE(numout,*) 'qsr : penetrative solar radiation forcing field written in ocean restart file ',   & 
     248            &                    'at it= ', kt,' date= ', ndastp 
     249         IF(lwp) WRITE(numout,*) '~~~~' 
     250         CALL iom_rstput( kt, nitrst, numrow, 'qsr_trd_hc_b', qsr_trd_hc_n ) 
     251         ! 
    227252      ENDIF 
    228253 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90

    r1975 r2068  
    2020   USE trdmod          ! ocean trends  
    2121   USE trdmod_oce      ! ocean variables trends 
     22   USE iom 
    2223   USE in_out_manager  ! I/O manager 
     24   USE restart         ! ocean restart 
    2325   USE prtctl          ! Print control 
    2426 
     
    132134         ENDIF 
    133135      ENDIF 
    134  
    135       !                             ! ---------------------- ! 
    136       IF( lk_vvl ) THEN             !  Variable Volume case  ! 
    137          !                          ! ---------------------- ! 
     136      !                                            ! ---------------------------------------- ! 
     137      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          ! 
     138         !                                         ! ---------------------------------------- ! 
     139         sbc_trd_hc_b(:,:) = sbc_trd_hc_n(:,:)                         ! Swap the ocean forcing fields except at nit000 
     140         IF ( .NOT. lk_vvl ) sbc_trd_sc_b(:,:)   = sbc_trd_sc_n(:,:) 
     141      ENDIF 
     142      !                                            ! ---------------------------------------- ! 
     143      IF( kt == nit000 ) THEN                      !   set the forcing field at nit000 - 1    ! 
     144         !                                         ! ---------------------------------------- ! 
     145         IF( ln_rstart .AND.    &                               !* Restart: read in restart file 
     146            & iom_varid( numror, 'sbc_trd_hc_b', ldstop = .FALSE. ) > 0 ) THEN  
     147            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file' 
     148            CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_hc_b', sbc_trd_hc_b )   ! before heat content sbc trend 
     149            CALL iom_get( numror, jpdom_autoglo, 'qsr_trd_hc_b', qsr_trd_hc_b )   ! before heat content trend due to Qsr flux 
     150            IF ( .NOT. lk_vvl ) THEN 
     151               CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_sc_b', sbc_trd_sc_b )   ! before salt content sbc trend 
     152            ENDIF 
     153         ENDIF 
     154      ENDIF 
     155      !                                            ! ---------------------- ! 
     156      IF( lk_vvl ) THEN                            !  Variable Volume case  ! 
     157         !                                         ! ---------------------- ! 
    138158!!gm BUG : in key_vvl emps must be modified to only include the salt flux due to sea-ice freezing/melting 
    139159!!gm       otherwise this flux will be missing  ==> modification required in limsbc,  limsbc_2 and CICE interface.s 
     
    161181            END DO 
    162182         ENDIF 
    163          !                          ! ---------------------- ! 
    164       ELSE                          !  Constant Volume case  ! 
    165          !                          ! ---------------------- ! 
     183         !                                         ! ---------------------- ! 
     184      ELSE                                         !  Constant Volume case  ! 
     185         !                                         ! ---------------------- ! 
    166186         IF ( neuler == 0 .AND. kt == nit000 ) THEN 
    167187            DO jj = 2, jpj 
     
    197217      ENDIF 
    198218 
     219      !                                            ! ---------------------------------------- ! 
     220      IF( lrst_oce ) THEN                          !      Write in the ocean restart file     ! 
     221         !                                         ! ---------------------------------------- ! 
     222         IF(lwp) WRITE(numout,*) 
     223         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   & 
     224            &                    'at it= ', kt,' date= ', ndastp 
     225         IF(lwp) WRITE(numout,*) '~~~~' 
     226         CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_hc_b', sbc_trd_hc_n ) 
     227         IF ( .NOT. lk_vvl ) CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_sc_b', sbc_trd_sc_n ) 
     228         ! 
     229      ENDIF 
     230 
    199231      IF( l_trdtra ) THEN           ! save the sbc trends for diagnostic 
    200232         ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:) 
     
    205237      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' sbc  - Ta: ', mask1=tmask,   & 
    206238         &                       tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    207       ! 
    208239   END SUBROUTINE tra_sbc 
    209240 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/istate.F90

    r1566 r2068  
    6666      !!---------------------------------------------------------------------- 
    6767      USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine) 
     68      ! - ML - needed for initialization of e3t_b 
     69      INTEGER  ::  jk     ! dummy loop indice 
    6870 
    6971      IF(lwp) WRITE(numout,*) 
     
    128130            &                                  gtu, gsu, gru, &  ! of t, s, rd at the bottom ocean level 
    129131            &                                  gtv, gsv, grv ) 
    130           
     132       
     133         ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here 
     134         IF( lk_vvl ) THEN 
     135            DO jk = 1, jpk 
     136               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 
     137            ENDDO 
     138         ENDIF 
    131139      ENDIF 
    132140      ! 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/step.F90

    r1975 r2068  
    195195      !  Ocean dynamics : ssh, wn, hdiv, rot                                 ! 
    196196      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     197 
    197198                         CALL ssh_wzv( kstp )         ! after ssh & vertical velocity 
    198199 
Note: See TracChangeset for help on using the changeset viewer.