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

Changeset 1870


Ignore:
Timestamp:
2010-05-12T17:36:00+02:00 (14 years ago)
Author:
gm
Message:

ticket: #663 step-1 : introduce the modified forcing term

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

Legend:

Unmodified
Added
Removed

  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynzdf_exp.F90

    r1537 r1870  
    44   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend 
    55   !!============================================================================== 
    6    !! History :      !  90-10  (B. Blanke)  Original code 
    7    !!                !  97-05  (G. Madec)  vertical component of isopycnal 
    8    !!           8.5  !  02-08  (G. Madec)  F90: Free form and module 
     6   !! History :  OPA  !  1990-10  (B. Blanke)  Original code 
     7   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal 
     8   !!   NEMO     1.0  !  1002-08  (G. Madec)  F90: Free form and module 
     9   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
    910   !!---------------------------------------------------------------------- 
    1011 
    … …  
    1314   !!                  sion using an explicit time-stepping scheme. 
    1415   !!---------------------------------------------------------------------- 
    15    !! * Modules used 
    1616   USE oce             ! ocean dynamics and tracers 
    1717   USE dom_oce         ! ocean space and time domain 
    … …  
    2424   PRIVATE 
    2525 
    26    !! * Routine accessibility 
    27    PUBLIC dyn_zdf_exp    ! called by step.F90 
     26   PUBLIC   dyn_zdf_exp   ! called by step.F90 
    2827 
    2928   !! * Substitutions 
    … …  
    3130#  include "vectopt_loop_substitute.h90" 
    3231   !!---------------------------------------------------------------------- 
    33    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     32   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    3433   !! $Id$ 
    3534   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    4746      !!      technique. The vertical diffusion of momentum is given by: 
    4847      !!         diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) ) 
    49       !!      Surface boundary conditions: wind stress input 
     48      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) 
    5049      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F90) 
    5150      !!      Add this trend to the general trend ua : 
    … …  
    5453      !! ** Action : - Update (ua,va) with the vertical diffusive trend 
    5554      !!--------------------------------------------------------------------- 
    56       !! * Arguments 
    57       INTEGER , INTENT( in ) ::   kt                           ! ocean time-step index 
    58       REAL(wp), INTENT( in ) ::   p2dt                         ! time-step  
    59  
    60       !! * Local declarations 
    61       INTEGER ::   ji, jj, jk, jl                              ! dummy loop indices 
    62       REAL(wp) ::   zrau0r, zlavmr, zua, zva                   ! temporary scalars 
    63       REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww     ! temporary workspace arrays 
     55      INTEGER , INTENT(in) ::   kt     ! ocean time-step index 
     56      REAL(wp), INTENT(in) ::   p2dt   ! time-step  
     57      !! 
     58      INTEGER ::   ji, jj, jk, jl                            ! dummy loop indices 
     59      REAL(wp) ::   zrau0r, zlavmr, zua, zva                 ! temporary scalars 
     60      REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww   ! 2D workspace 
    6461      !!---------------------------------------------------------------------- 
    6562 
    6663      IF( kt == nit000 ) THEN 
    6764         IF(lwp) WRITE(numout,*) 
    68          IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion explicit operator' 
     65         IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion - explicit operator' 
    6966         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
    7067      ENDIF 
    7168 
    72       ! Local constant initialization 
    73       ! ----------------------------- 
    74       zrau0r = 1. / rau0                                   ! inverse of the reference density 
    75       zlavmr = 1. / float( nn_zdfexp )                      ! inverse of the number of sub time step 
     69      zrau0r = 1. / rau0                              ! Local constant initialization 
     70      zlavmr = 1. / REAL( nn_zdfexp ) 
    7671 
    77       !                                                ! =============== 
    78       DO jj = 2, jpjm1                                 !  Vertical slab 
    79          !                                             ! =============== 
    80  
    81          ! Surface boundary condition 
    82          DO ji = 2, jpim1 
    83             zwy(ji,1) = utau(ji,jj) * zrau0r 
    84             zww(ji,1) = vtau(ji,jj) * zrau0r 
     72      !                                               ! =============== 
     73      DO jj = 2, jpjm1                                !  Vertical slab 
     74         !                                            ! =============== 
     75         DO ji = 2, jpim1         ! Surface boundary condition 
     76            zwy(ji,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * zrau0r 
     77            zww(ji,1) = ( vtau_b(ji,jj) + vtau(ji,jj) ) * zrau0r 
    8578         END DO   
    86  
    87          ! Initialization of x, z and contingently trends array 
    88          DO jk = 1, jpk 
     79         DO jk = 1, jpk         ! Initialization of x, z and contingently trends array 
    8980            DO ji = 2, jpim1 
    9081               zwx(ji,jk) = ub(ji,jj,jk) 
    … …  
    9283            END DO   
    9384         END DO   
    94  
    95          ! Time splitting loop 
    96          DO jl = 1, nn_zdfexp 
    97  
    98             ! First vertical derivative 
    99             DO jk = 2, jpk 
     85         ! 
     86         DO jl = 1, nn_zdfexp     ! Time splitting loop 
     87            ! 
     88            DO jk = 2, jpk            ! First vertical derivative 
    10089               DO ji = 2, jpim1 
    10190                  zwy(ji,jk) = avmu(ji,jj,jk) * ( zwx(ji,jk-1) - zwx(ji,jk) ) / fse3uw(ji,jj,jk)  
    … …  
    10392               END DO   
    10493            END DO   
    105  
    106             ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp 
    107             DO jk = 1, jpkm1 
     94            DO jk = 1, jpkm1          ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp 
    10895               DO ji = 2, jpim1 
    109                   zua = zlavmr*( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk) 
    110                   zva = zlavmr*( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk) 
     96                  zua = zlavmr * ( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk) 
     97                  zva = zlavmr * ( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk) 
    11198                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    11299                  va(ji,jj,jk) = va(ji,jj,jk) + zva 
    113  
    114                   zwx(ji,jk) = zwx(ji,jk) + p2dt*zua*umask(ji,jj,jk) 
    115                   zwz(ji,jk) = zwz(ji,jk) + p2dt*zva*vmask(ji,jj,jk) 
     100                  ! 
     101                  zwx(ji,jk) = zwx(ji,jk) + p2dt * zua * umask(ji,jj,jk) 
     102                  zwz(ji,jk) = zwz(ji,jk) + p2dt * zva * vmask(ji,jj,jk) 
    116103               END DO   
    117104            END DO   
    118  
     105            ! 
    119106         END DO   
    120  
    121          !                                             ! =============== 
    122       END DO                                           !   End of slab 
    123       !                                                ! =============== 
    124  
     107         !                                            ! =============== 
     108      END DO                                          !   End of slab 
     109      !                                               ! =============== 
    125110   END SUBROUTINE dyn_zdf_exp 
    126111 

  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r1662 r1870  
    44   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend 
    55   !!============================================================================== 
    6    !! History :      !  90-10  (B. Blanke)  Original code 
    7    !!                !  97-05  (G. Madec)  vertical component of isopycnal 
    8    !!           8.5  !  02-08  (G. Madec)  F90: Free form and module 
     6   !! History :  OPA  !  1990-10  (B. Blanke)  Original code 
     7   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal 
     8   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module 
     9   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
    910   !!---------------------------------------------------------------------- 
    1011 
    … …  
    1314   !!                  sion using a implicit time-stepping. 
    1415   !!---------------------------------------------------------------------- 
    15    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
    16    !! $Id$ 
    17    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    18    !!---------------------------------------------------------------------- 
    19    !! * Modules used 
    2016   USE oce             ! ocean dynamics and tracers 
    2117   USE dom_oce         ! ocean space and time domain 
    … …  
    2824   PRIVATE 
    2925 
    30    !! * Routine accessibility 
    31    PUBLIC dyn_zdf_imp    ! called by step.F90 
     26   PUBLIC   dyn_zdf_imp   ! called by step.F90 
    3227 
    3328   !! * Substitutions 
    … …  
    3530#  include "vectopt_loop_substitute.h90" 
    3631   !!---------------------------------------------------------------------- 
    37    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     32   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    3833   !! $Id$ 
    39    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     34   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4035   !!---------------------------------------------------------------------- 
    4136 
    4237CONTAINS 
    43  
    4438 
    4539   SUBROUTINE dyn_zdf_imp( kt, p2dt ) 
    … …  
    5448      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) ) 
    5549      !!      backward time stepping 
    56       !!      Surface boundary conditions: wind stress input 
     50      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) 
    5751      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F) 
    5852      !!      Add this trend to the general trend ua : 
    5953      !!         ua = ua + dz( avmu dz(u) ) 
    6054      !! 
    61       !! ** Action : - Update (ua,va) arrays with the after vertical diffusive 
    62       !!               mixing trend. 
     55      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend. 
    6356      !!--------------------------------------------------------------------- 
    64       !! * Modules used 
    65       USE oce, ONLY :  zwd   => ta,   &                ! use ta as workspace 
    66                        zws   => sa                     ! use sa as workspace 
    67  
    68       !! * Arguments 
    69       INTEGER , INTENT( in ) ::   kt                   ! ocean time-step index 
    70       REAL(wp), INTENT( in ) ::  p2dt                  ! vertical profile of tracer time-step 
    71  
    72       !! * Local declarations 
    73       INTEGER ::   ji, jj, jk                          ! dummy loop indices 
    74       REAL(wp) ::   zrau0r, z2dtf, zcoef, zzws, zrhs   ! temporary scalars 
    75       REAL(wp) ::   zzwi                               ! temporary scalars 
    76       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays 
     57      USE oce, ONLY :  zwd   => ta      ! use ta as workspace 
     58      USE oce, ONLY :  zws   => sa      ! use sa as workspace 
     59      !! 
     60      INTEGER , INTENT( in ) ::   kt    ! ocean time-step index 
     61      REAL(wp), INTENT( in ) ::  p2dt   ! vertical profile of tracer time-step 
     62      !! 
     63      INTEGER  ::   ji, jj, jk             ! dummy loop indices 
     64      REAL(wp) ::   zrau0r, z2dtf, zcoef   ! temporary scalars 
     65      REAL(wp) ::   zzwi, zzws, zrhs       ! temporary scalars 
     66      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! 3D workspace 
    7767      !!---------------------------------------------------------------------- 
    7868 
    … …  
    9585      ! friction velocity in dyn_bfr using values from the previous timestep. There 
    9686      ! is no need to include these in the implicit calculation. 
    97       DO jk = 1, jpkm1 
     87      ! 
     88      DO jk = 1, jpkm1        ! Matrix 
    9889         DO jj = 2, jpjm1  
    9990            DO ji = fs_2, fs_jpim1   ! vector opt. 
    10091               zcoef = - p2dt / fse3u(ji,jj,jk) 
    101                zzwi          = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
    102                zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 
    103                zzws          = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
    104                zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 
     92               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
     93               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk) 
     94               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
     95               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1) 
    10596               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws 
    10697            END DO 
    10798         END DO 
    10899      END DO 
    109  
    110       ! Surface boudary conditions 
    111       DO jj = 2, jpjm1  
     100      DO jj = 2, jpjm1        ! Surface boudary conditions 
    112101         DO ji = fs_2, fs_jpim1   ! vector opt. 
    113102            zwi(ji,jj,1) = 0. 
    … …  
    130119      !   The solution (the after velocity) is in ua 
    131120      !----------------------------------------------------------------------- 
    132  
    133       ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k) 
    134       DO jk = 2, jpkm1 
     121      ! 
     122      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    135123         DO jj = 2, jpjm1    
    136124            DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    139127         END DO 
    140128      END DO 
    141  
    142       ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1 
    143       DO jj = 2, jpjm1    
    144          DO ji = fs_2, fs_jpim1   ! vector opt. 
    145 !!! change les resultats (derniers digit, pas significativement + rapide 1* de moins) 
    146 !!!         ua(ji,jj,1) = ub(ji,jj,1)  & 
    147 !!!                      + p2dt * ( ua(ji,jj,1) + utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) ) 
    148             z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 ) 
    149             ua(ji,jj,1) = ub(ji,jj,1)  & 
    150                          + p2dt *  ua(ji,jj,1) + z2dtf * utau(ji,jj) 
     129      ! 
     130      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
     131         DO ji = fs_2, fs_jpim1   ! vector opt. 
     132            ua(ji,jj,1) = ub(ji,jj,1)   & 
     133               &        + p2dt * (  ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 )  ) 
    151134         END DO 
    152135      END DO 
    … …  
    159142         END DO 
    160143      END DO 
    161  
    162       ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 
    163       DO jj = 2, jpjm1    
     144      ! 
     145      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  == 
    164146         DO ji = fs_2, fs_jpim1   ! vector opt. 
    165147            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    … …  
    173155         END DO 
    174156      END DO 
    175  
    176       ! Normalization to obtain the general momentum trend ua 
    177       DO jk = 1, jpkm1 
     157      ! 
     158      DO jk = 1, jpkm1        !==  Normalization to obtain the general momentum trend ua  == 
    178159         DO jj = 2, jpjm1    
    179160            DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    192173      ! friction velocity in dyn_bfr using values from the previous timestep. There 
    193174      ! is no need to include these in the implicit calculation. 
    194       DO jk = 1, jpkm1 
     175      ! 
     176      DO jk = 1, jpkm1        ! Matrix 
    195177         DO jj = 2, jpjm1    
    196178            DO ji = fs_2, fs_jpim1   ! vector opt. 
    197179               zcoef = -p2dt / fse3v(ji,jj,jk) 
    198                zzwi          = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
     180               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
    199181               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk) 
    200                zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
     182               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
    201183               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1) 
    202184               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws 
    … …  
    204186         END DO 
    205187      END DO 
    206  
    207       ! Surface boudary conditions 
    208       DO jj = 2, jpjm1    
     188      DO jj = 2, jpjm1        ! Surface boudary conditions 
    209189         DO ji = fs_2, fs_jpim1   ! vector opt. 
    210190            zwi(ji,jj,1) = 0.e0 
    … …  
    223203      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
    224204      ! 
    225       !   m is decomposed in the product of an upper and lower triangular 
    226       !   matrix 
     205      !   m is decomposed in the product of an upper and lower triangular matrix 
    227206      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    228207      !   The solution (after velocity) is in 2d array va 
    229208      !----------------------------------------------------------------------- 
    230  
    231       ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k) 
    232       DO jk = 2, jpkm1 
     209      ! 
     210      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    233211         DO jj = 2, jpjm1    
    234212            DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    237215         END DO 
    238216      END DO 
    239  
    240       ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1 
    241       DO jj = 2, jpjm1 
    242          DO ji = fs_2, fs_jpim1   ! vector opt. 
    243 !!! change les resultats (derniers digit, pas significativement + rapide 1* de moins) 
    244 !!!         va(ji,jj,1) = vb(ji,jj,1)  & 
    245 !!!                      + p2dt * ( va(ji,jj,1) + vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) ) 
    246             z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 ) 
    247             va(ji,jj,1) = vb(ji,jj,1)  & 
    248                          + p2dt * va(ji,jj,1) + z2dtf * vtau(ji,jj) 
     217      ! 
     218      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
     219         DO ji = fs_2, fs_jpim1   ! vector opt. 
     220            va(ji,jj,1) = vb(ji,jj,1)   & 
     221               &        + p2dt * (  va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 )  ) 
    249222         END DO 
    250223      END DO 
    … …  
    257230         END DO 
    258231      END DO 
    259  
    260       ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 
    261       DO jj = 2, jpjm1    
     232      ! 
     233      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  == 
    262234         DO ji = fs_2, fs_jpim1   ! vector opt. 
    263235            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    … …  
    271243         END DO 
    272244      END DO 
    273  
    274       ! Normalization to obtain the general momentum trend va 
    275       DO jk = 1, jpkm1 
     245      ! 
     246      DO jk = 1, jpkm1     !==  Normalization to obtain the general momentum trend va  == 
    276247         DO jj = 2, jpjm1    
    277248            DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    280251         END DO 
    281252      END DO 
    282  
     253      ! 
    283254   END SUBROUTINE dyn_zdf_imp 
    284255 

  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90

    r1792 r1870  
    55   !!============================================================================== 
    66   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code 
     7   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA  
    78   !!---------------------------------------------------------------------- 
    89 
    … …  
    3738#  include "domzgr_substitute.h90" 
    3839#  include "vectopt_loop_substitute.h90" 
    39  
    40    !!---------------------------------------------------------------------- 
    41    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     40   !!---------------------------------------------------------------------- 
     41   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    4242   !! $Id$ 
    4343   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    5353      !!              and update the now vertical coordinate (lk_vvl=T). 
    5454      !! 
    55       !! ** Method  : - 
    56       !! 
    57       !!              - Using the incompressibility hypothesis, the vertical  
     55      !! ** Method  : - Using the incompressibility hypothesis, the vertical  
    5856      !!      velocity is computed by integrating the horizontal divergence   
    5957      !!      from the bottom to the surface minus the scale factor evolution. 
    … …  
    6260      !! ** action  :   ssha    : after sea surface height 
    6361      !!                wn      : now vertical velocity 
    64       !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height 
    65       !!                          at u-, v-, f-point s 
    66       !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 
     62      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T) 
     63      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points 
     64      !! 
     65      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    6766      !!---------------------------------------------------------------------- 
    6867      USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace 
    … …  
    7877 
    7978      IF( kt == nit000 ) THEN 
     79         ! 
    8080         IF(lwp) WRITE(numout,*) 
    8181         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' 
    … …  
    111111      ENDIF 
    112112 
    113       !                                           !------------------------------! 
    114       IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case) 
    115       !                                           !------------------------------! 
     113      !                                           !------------------------------------------! 
     114      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case) 
     115         !                                        !------------------------------------------! 
    116116         DO jk = 1, jpkm1 
    117             fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays 
     117            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays 
    118118            fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 
    119119            fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 
    120120            ! 
    121             fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays 
     121            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays 
    122122            fse3u (:,:,jk) = fse3u_n (:,:,jk) 
    123123            fse3v (:,:,jk) = fse3v_n (:,:,jk) 
    … …  
    127127            fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 
    128128         END DO 
    129          !                                             ! now ocean depth (at u- and v-points) 
    130          hu(:,:) = hu_0(:,:) + sshu_n(:,:) 
     129         ! 
     130         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points) 
    131131         hv(:,:) = hv_0(:,:) + sshv_n(:,:) 
    132          !                                             ! now masked inverse of the ocean depth (at u- and v-points) 
     132         !                                            ! now masked inverse of the ocean depth (at u- and v-points) 
    133133         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) ) 
    134134         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) ) 
    … …  
    136136      ENDIF 
    137137 
    138                          CALL div_cur( kt )            ! Horizontal divergence & Relative vorticity 
    139       IF( n_cla == 1 )   CALL div_cla( kt )            ! Cross Land Advection (Update Hor. divergence) 
    140  
    141       ! set time step size (Euler/Leapfrog) 
    142       z2dt = 2. * rdt  
     138                         CALL div_cur( kt )           ! Horizontal divergence & Relative vorticity 
     139      IF( n_cla == 1 )   CALL div_cla( kt )           ! Cross Land Advection (Update Hor. divergence) 
     140 
     141      z2dt = 2. * rdt                                 ! set time step size (Euler/Leapfrog) 
    143142      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt 
    144143 
    145       zraur = 1. / rau0 
    146  
    147144      !                                           !------------------------------! 
    148145      !                                           !   After Sea Surface Height   ! 
    149146      !                                           !------------------------------! 
     147      zraur = 0.5 / rau0 
     148      ! 
    150149      zhdiv(:,:) = 0.e0 
    151150      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
    152151        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 
    153152      END DO 
    154  
    155153      !                                                ! Sea surface elevation time stepping 
    156       ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     154      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    157155 
    158156#if defined key_obc 
    159       IF ( Agrif_Root() ) THEN  
     157      IF( Agrif_Root() ) THEN  
    160158         ssha(:,:) = ssha(:,:) * obctmsk(:,:) 
    161          CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm) 
     159         CALL lbc_lnk( ssha, 'T', 1. )                ! absolutly compulsory !! (jmm) 
    162160      ENDIF 
    163161#endif 
    164  
    165162      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only) 
    166163      IF( lk_vvl ) THEN                                ! (required only in key_vvl case) 
    … …  
    186183      !                                           !     Now Vertical Velocity    ! 
    187184      !                                           !------------------------------! 
    188       !                                                ! integrate from the bottom the hor. divergence 
    189       DO jk = jpkm1, 1, -1 
    190          wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   & 
    191               &                    - (  fse3t_a(:,:,jk)                   & 
    192               &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 
     185      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
     186         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)     & 
     187            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 
    193188      END DO 
    194       ! 
     189 
     190      !                                           !------------------------------! 
     191      !                                           !           outputs            ! 
     192      !                                           !------------------------------! 
    195193      CALL iom_put( "woce", wn                    )   ! vertical velocity 
    196194      CALL iom_put( "ssh" , sshn                  )   ! sea surface height 
    197195      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height 
    198       IF( lk_diaar5 ) THEN 
     196      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value 
    199197         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) 
    200198         DO jk = 1, jpk 
    201199            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 
    202200         END DO 
    203          CALL iom_put( "w_masstr" , z3d                     )   !           vertical mass transport 
    204          CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )   ! square of vertical mass transport 
    205       ENDIF 
     201         CALL iom_put( "w_masstr" , z3d                     )   
     202         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
     203      ENDIF 
     204      ! 
     205      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 ) 
    206206      ! 
    207207   END SUBROUTINE ssh_wzv 
    … …  
    216216      !!              ssha  already computed in ssh_wzv   
    217217      !! 
    218       !! ** Method  : - apply Asselin time fiter to now ssh and swap : 
    219       !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 
    220       !!             sshn = ssha 
     218      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing 
     219      !!              from the filter, see Leclair and Madec 2010) and swap : 
     220      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 
     221      !!                            - atfp * rdt * ( emp_b - emp ) / rau0 
     222      !!                sshn = ssha 
    221223      !! 
    222224      !! ** action  : - sshb, sshn   : before & now sea surface height 
    223225      !!                               ready for the next time step 
    224       !!---------------------------------------------------------------------- 
    225       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    226       !! 
    227       INTEGER  ::   ji, jj               ! dummy loop indices 
     226      !! 
     227      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     228      !!---------------------------------------------------------------------- 
     229      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     230      !! 
     231      INTEGER  ::   ji, jj   ! dummy loop indices 
     232      REAL(wp) ::   zec      ! temporary scalare 
    228233      !!---------------------------------------------------------------------- 
    229234 
    … …  
    234239      ENDIF 
    235240 
    236       ! Time filter and swap of the ssh 
    237       ! ------------------------------- 
    238  
    239       IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points 
    240          !                   ! ---------------------- ! 
     241      !                       !--------------------------! 
     242      IF( lk_vvl ) THEN       !  Variable volume levels  !   ssh at t-, u-, v, f-points 
     243         !                    !--------------------------! 
    241244         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter 
    242245            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now) 
    … …  
    247250            DO jj = 1, jpj 
    248251               DO ji = 1, jpi                                ! before <-- now filtered 
    249                   sshb  (ji,jj) = sshn(ji,jj)  + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) ) 
     252                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) ) 
    250253                  sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 
    251254                  sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) 
    … …  
    258261            END DO 
    259262         ENDIF 
    260          ! 
    261       ELSE                   ! fixed levels :   ssh at t-point only 
    262          !                   ! ------------ ! 
     263         !                    !--------------------------! 
     264      ELSE                    !        fixed levels      !   ssh at t-point only 
     265         !                    !--------------------------! 
    263266         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter 
    264267            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now) 
    265268            ! 
    266269         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap 
     270            zec = atfp * rdt / rau0 
    267271            DO jj = 1, jpj 
    268272               DO ji = 1, jpi                                ! before <-- now filtered 
    269                   sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )     
     273                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   & 
     274                     &                      - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 
    270275                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after 
    271276               END DO 
    … …  
    274279         ! 
    275280      ENDIF 
    276       ! 
    277       IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 ) 
     281 
     282      ! time filter with global conservation correction and array swap 
     283      sshbb(:,:) = sshb(:,:) 
     284      sshb (:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + zssha(:,:) )    & 
     285           &       - zfact  *  
     286      sshn (:,:) = zssha(:,:) 
     287      empb (:,:) = emp(:,:) 
     288 
     289 
     290      ! 
     291      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 ) 
    278292      ! 
    279293   END SUBROUTINE ssh_nxt 

  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r1705 r1870  
    44   !! Surface module :   variables defined in core memory  
    55   !!====================================================================== 
    6    !! History :  3.0   !  2006-06  (G. Madec)  Original code 
    7    !!             -    !  2008-08  (G. Madec)  namsbc moved from sbcmod 
     6   !! History :  3.0  !  2006-06  (G. Madec)  Original code 
     7   !!             -   !  2008-08  (G. Madec)  namsbc moved from sbcmod 
     8   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
    89   !!---------------------------------------------------------------------- 
    910   USE par_oce          ! ocean parameters 
    … …  
    2526   LOGICAL , PUBLIC ::   ln_ssr      = .FALSE.   !: Sea Surface restoring on SST and/or SSS       
    2627   INTEGER , PUBLIC ::   nn_ice      = 0         !: flag on ice in the surface boundary condition (=0/1/2/3) 
    27    INTEGER , PUBLIC ::   nn_fwb      = 0         !: FreshWater Budget:  
     28   INTEGER , PUBLIC ::   nn_fwb      = 0         !: FreshWater Budget control :  
    2829   !                                             !:  = 0 unchecked  
    2930   !                                             !:  = 1 global mean of e-p-r set to zero at each nn_fsbc time step 
    3031   !                                             !:  = 2 annual global mean of e-p-r set to zero 
    31    INTEGER , PUBLIC ::   nn_ico_cpl  = 0          !: ice-ocean coupling indicator 
     32   INTEGER , PUBLIC ::   nn_ico_cpl  = 0         !: ice-ocean coupling indicator 
    3233   !                                             !:  = 0   LIM-3 old case 
    3334   !                                             !:  = 1   stresses computed using now ocean velocity 
    … …  
    3738   !!              Ocean Surface Boundary Condition fields 
    3839   !!---------------------------------------------------------------------- 
    39    LOGICAL , PUBLIC ::   lhftau = .FALSE.              !: HF tau contribution: mean of stress module - module of the mean stress 
    40    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau      !: sea surface i-stress (ocean referential)     [N/m2] 
    41    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau      !: sea surface j-stress (ocean referential)     [N/m2] 
    42    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   taum      !: module of sea surface stress (at T-point)    [N/m2]  
    43    !! wndm is used only in PISCES to compute gases exchanges at the surface of the free ocean or in the leads in sea-ice parts 
    44    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   wndm      !: wind speed module at T-point (=|U10m-Uoce|)  [m/s]  
    45    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr       !: sea heat flux:     solar                     [W/m2] 
    46    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns       !: sea heat flux: non solar                     [W/m2] 
    47    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr_tot   !: total     solar heat flux (over sea and ice) [W/m2] 
    48    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns_tot   !: total non solar heat flux (over sea and ice) [W/m2] 
    49    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp       !: freshwater budget: volume flux               [Kg/m2/s] 
    50    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps      !: freshwater budget: concentration/dillution   [Kg/m2/s] 
    51    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot   !: total evaporation - (liquid + solid) precpitation over oce and ice 
    52    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tprecip   !: total precipitation           [Kg/m2/s] 
    53    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sprecip   !: solid precipitation           [Kg/m2/s] 
    54 !!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rrunoff       !: runoff 
    55 !!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   calving       !: calving 
    56    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr_i      !: ice fraction  (between 0 to 1)               - 
     40   LOGICAL , PUBLIC ::   lhftau = .FALSE.              !: HF tau used in TKE: mean(stress module) - module(mean stress) 
     41   !!                                   !!   now    ! before   !! 
     42   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2] 
     43   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2] 
     44   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   taum              !: module of sea surface stress (at T-point)    [N/m2]  
     45   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s] 
     46   !! wndm is used only in PISCES to compute surface gases exchanges in ice-free ocean or leads 
     47   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr    , qsr_b    !: sea heat flux:     solar                     [W/m2] 
     48   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns    , qns_b    !: sea heat flux: non solar                     [W/m2] 
     49   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr_tot           !: total     solar heat flux (over sea and ice) [W/m2] 
     50   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns_tot           !: total non solar heat flux (over sea and ice) [W/m2] 
     51   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp    , emp_b    !: freshwater budget: volume flux               [Kg/m2/s] 
     52   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps   , emps_b   !: freshwater budget: concentration/dillution   [Kg/m2/s] 
     53   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot           !: total E-P-R over ocean and ice               [Kg/m2/s] 
     54   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tprecip           !: total precipitation                          [Kg/m2/s] 
     55   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sprecip           !: solid precipitation                          [Kg/m2/s] 
     56   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr_i              !: ice fraction = 1 - lead fraction      (between 0 to 1) 
    5757#if defined key_cpl_carbon_cycle 
    58    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   atm_co2   !: atmospheric pCO2                             [ppm] 
     58   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   atm_co2           !: atmospheric pCO2                             [ppm] 
    5959#endif 
     60!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rrunoff        !: runoff 
     61!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   calving        !: calving 
    6062 
    6163   !!---------------------------------------------------------------------- 
    … …  
    7072 
    7173   !!---------------------------------------------------------------------- 
    72    !!   OPA 9.0 , LOCEAN-IPSL (2006)  
    73    !! $ Id: $ 
     74   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
     75   !! $Id$ 
    7476   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    7577   !!====================================================================== 
    76  
    7778END MODULE sbc_oce 

  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90

    r1792 r1870  
    44   !! Surface module :  provide to the ocean its surface boundary condition 
    55   !!====================================================================== 
    6    !! History :  3.0   !  07-2006  (G. Madec)  Original code 
    7    !!             -    !  08-2008  (S. Masson, E. .... ) coupled interface 
     6   !! History :  3.0  !  2006-07  (G. Madec)  Original code 
     7   !!            3.1  !  2008-08  (S. Masson, E. Maisonnave, G. Madec) coupled interface 
     8   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
    89   !!---------------------------------------------------------------------- 
    910 
    … …  
    4950#  include "domzgr_substitute.h90" 
    5051   !!---------------------------------------------------------------------- 
    51    !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     52   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    5253   !! $Id$ 
    5354   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    8687!!gm  IF( lk_sbc_cpl       ) THEN   ;   ln_cpl      = .TRUE.   ;   ELSE   ;   ln_cpl      = .FALSE.   ;   ENDIF 
    8788 
    88       IF ( Agrif_Root() ) THEN 
     89      IF( Agrif_Root() ) THEN 
    8990        IF( lk_lim2 )            nn_ice      = 2 
    9091        IF( lk_lim3 )            nn_ice      = 3 
    … …  
    179180      !!                CAUTION : never mask the surface stress field (tke sbc) 
    180181      !! 
    181       !! ** Action  : - set the ocean surface boundary condition, i.e.   
    182       !!                utau, vtau, qns, qsr, emp, emps, qrp, erp 
     182      !! ** Action  : - set the ocean surface boundary condition at before and now  
     183      !!                time step, i.e.   
     184      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, emps_b, qrp_b, erp_b 
     185      !!                utau  , vtau  , qns  , qsr  , emp  , emps  , qrp  , erp 
    183186      !!              - updte the ice fraction : fr_i 
    184187      !!---------------------------------------------------------------------- 
    … …  
    186189      !!--------------------------------------------------------------------- 
    187190 
    188       CALL iom_setkt( kt + nn_fsbc - 1 )         !  in sbc, iom_put is called every nn_fsbc time step 
    189       ! 
    190       ! ocean to sbc mean sea surface variables (ss._m) 
    191       ! --------------------------------------- 
    192       CALL sbc_ssm( kt )                         ! sea surface mean currents (at U- and V-points),  
    193       !                                          ! temperature and salinity (at T-point) over nf_sbc time-step 
    194       !                                          ! (i.e. sst_m, sss_m, ssu_m, ssv_m) 
    195  
    196       ! sbc formulation 
    197       ! --------------- 
    198           
    199       SELECT CASE( nsbc )                        ! Compute ocean surface boundary condition 
    200       !                                          ! (i.e. utau,vtau, qns, qsr, emp, emps) 
     191      !                                            ! ---------------------------------------- ! 
     192      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          ! 
     193         !                                         ! ---------------------------------------- ! 
     194         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields 
     195         utau_b(:,:) = utau(:,:)                         ! (except at nitOOO where before fields 
     196         qns_b (:,:) = qns (:,:)                         !  are set the end of the routine) 
     197         qsr_b (:,:) = qsr (:,:) 
     198         emp_b (:,:) = emp (:,:) 
     199         emps_b(:,:) = emps(:,:) 
     200      ENDIF 
     201 
     202      !                                            ! ---------------------------------------- ! 
     203      !                                            !        forcing field computation         ! 
     204      !                                            ! ---------------------------------------- ! 
     205 
     206      CALL iom_setkt( kt + nn_fsbc - 1 )                 ! in sbc, iom_put is called every nn_fsbc time step 
     207      ! 
     208      CALL sbc_ssm( kt )                                 ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) 
     209      !                                                  ! averaged over nf_sbc time-step 
     210 
     211                                                   !==  sbc formulation  ==! 
     212                                                             
     213      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition 
     214      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, emps) 
    201215      CASE(  0 )   ;   CALL sbc_gyre    ( kt )                    ! analytical formulation : GYRE configuration 
    202216      CASE(  1 )   ;   CALL sbc_ana     ( kt )                    ! analytical formulation : uniform sbc 
    … …  
    214228      END SELECT 
    215229 
    216       ! Misc. Options 
    217       ! ------------- 
     230      !                                            !==  Misc. Options  ==! 
    218231 
    219232!!gm  IF( ln_dm2dc       )   CALL sbc_dcy( kt )                 ! Daily mean qsr distributed over the Diurnal Cycle 
    … …  
    236249      !                                                         ! (update freshwater fluxes) 
    237250      ! 
     251      !                                                ! ---------------------------------------- ! 
     252      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    ! 
     253         !                                             ! ---------------------------------------- ! 
     254         IF( ln_rstart .AND.    &                               !* Restart: read in restart file 
     255            & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN  
     256            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file' 
     257            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before i-stress  (U-point) 
     258            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before j-stress  (V-point) 
     259            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point) 
     260            CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  )   ! before     solar heat flux (T-point) 
     261            CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b  )   ! before     freshwater flux (T-point) 
     262            CALL iom_get( numror, jpdom_autoglo, 'emps_b', emp_b  )   ! before C/D freshwater flux (T-point) 
     263            ! 
     264         ELSE                                                   !* no restart: set from nit000 values 
     265            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000' 
     266            utau_b(:,:) = utau(:,:)  
     267            utau_b(:,:) = utau(:,:) 
     268            qns_b (:,:) = qns (:,:) 
     269            qsr_b (:,:) = qsr (:,:) 
     270            emp_b (:,:) = emp (:,:) 
     271            emps_b(:,:) = emps(:,:) 
     272         ENDIF 
     273      ENDIF 
     274 
     275      !                                                ! ---------------------------------------- ! 
     276      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     ! 
     277         !                                             ! ---------------------------------------- ! 
     278         IF(lwp) WRITE(numout,*) 
     279         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   & 
     280            &                    'at it= ', kt,' date= ', ndastp 
     281         IF(lwp) WRITE(numout,*) '~~~~' 
     282         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )    !  
     283         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , vtau ) 
     284         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  ) 
     285         CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  ) 
     286         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  ) 
     287         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emp  ) 
     288         ! 
     289      ENDIF 
     290 
     291      !                                                ! ---------------------------------------- ! 
     292      !                                                !        Outputs and control print         ! 
     293      !                                                ! ---------------------------------------- ! 
    238294      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    239295         CALL iom_put( "emp"    , emp       )                   ! upward water flux 
    … …  
    242298         CALL iom_put( "qns"    , qns       )                   ! solar heat flux    moved after the call to iom_setkt) 
    243299         CALL iom_put( "qsr"    ,       qsr )                   ! solar heat flux    moved after the call to iom_setkt) 
    244          IF(  nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )  ! ice fraction  
     300         IF( nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction  
    245301      ENDIF 
    246302      ! 

  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90

    r1847 r1870  
    1515   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf 
    1616   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option 
    17    !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter 
     17   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA 
    1818   !!---------------------------------------------------------------------- 
    1919 
    … …  
    8181      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
    8282      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
     83      !! 
     84      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    8385      !!---------------------------------------------------------------------- 
    8486      USE oce, ONLY :    ztrdt => ua   ! use ua as 3D workspace    
    … …  
    9193      !!---------------------------------------------------------------------- 
    9294 
    93       IF( kt == nit000 ) THEN 
     95      IF( kt == nit000 ) THEN          !==  initialisation  ==! 
    9496         IF(lwp) WRITE(numout,*) 
    9597         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 
    … …  
    99101         r1_2bcp = 1.e0 - 2.e0 * rbcp 
    100102      ENDIF 
    101  
    102       ! Update after tracer on domain lateral boundaries 
     103      !                                      ! set time step size (Euler/Leapfrog) 
     104      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler) 
     105      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog) 
     106      ENDIF 
     107 
     108      !                                !==  Update after tracer on domain lateral boundaries  ==! 
    103109      !  
    104       CALL lbc_lnk( ta, 'T', 1. )      ! local domain boundaries  (T-point, unchanged sign) 
     110      CALL lbc_lnk( ta, 'T', 1. )            ! local domain boundaries  (T-point, unchanged sign) 
    105111      CALL lbc_lnk( sa, 'T', 1. ) 
    106112      ! 
    107113#if defined key_obc 
    108       CALL obc_tra( kt )               ! OBC open boundaries 
     114      CALL obc_tra( kt )                     ! OBC open boundaries 
    109115#endif 
    110116#if defined key_bdy 
    111       CALL bdy_tra( kt )               ! BDY open boundaries 
     117      CALL bdy_tra( kt )                     ! BDY open boundaries 
    112118#endif 
    113119#if defined key_agrif 
    114       CALL Agrif_tra                   ! AGRIF zoom boundaries 
     120      CALL Agrif_tra                         ! AGRIF zoom boundaries 
    115121#endif 
    116122  
    117       ! set time step size (Euler/Leapfrog) 
    118       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler) 
    119       ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog) 
    120       ENDIF 
    121  
    122       ! trends computation initialisation 
    123       IF( l_trdtra ) THEN              ! store now fields before applying the Asselin filter 
     123      IF( l_trdtra ) THEN              ! trends computation: store now fields before applying the Asselin filter 
    124124         ztrdt(:,:,:) = tn(:,:,:)       
    125125         ztrds(:,:,:) = sn(:,:,:) 
    126126      ENDIF 
    127127 
    128       ! Leap-Frog + Asselin filter time stepping 
     128      !                                !==  modifed Leap-Frog + Asselin filter (modified LF-RA) scheme  ==! 
    129129      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt )      ! variable volume level (vvl) 
    130130      ELSE                  ;   CALL tra_nxt_fix( kt )      ! fixed    volume level 
    … …  
    132132 
    133133#if defined key_agrif 
    134       ! Update tracer at AGRIF zoom boundaries 
    135       IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only 
     134      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )   ! Update tracer at AGRIF zoom boundaries (children only) 
    136135#endif       
    137136 
    138       ! trends computation 
    139       IF( l_trdtra ) THEN              ! trend of the Asselin filter (tb filtered - tb)/dt      
     137      IF( l_trdtra ) THEN              ! trends computation: trend of the Asselin filter (tb filtered - tb)/dt      
    140138         DO jk = 1, jpkm1 
    141139            zfact = 1.e0 / r2dt_t(jk)              
    … …  
    176174      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
    177175      !!---------------------------------------------------------------------- 
    178       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    179       !! 
    180       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    181       REAL(wp) ::   ztm, ztf     ! temporary scalars 
    182       REAL(wp) ::   zsm, zsf     !    -         - 
     176      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     177      !! 
     178      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
     179      REAL(wp) ::   ztm, ztf, zac   ! temporary scalars 
     180      REAL(wp) ::   zsm, zsf        !    -         - 
    183181      !!---------------------------------------------------------------------- 
    184182 
    … …  
    225223               sn(:,:,jk) = sa(:,:,jk) 
    226224            END DO 
    227          ELSE                                             ! general case (Leapfrog + Asselin filter 
     225         ELSE                                             ! general case (Leapfrog + Asselin filter) 
    228226            DO jk = 1, jpkm1 
    229227               DO jj = 1, jpj 
    … …  
    239237            END DO 
    240238         ENDIF 
     239      ENDIF 
     240      ! 
     241!!gm from Matthieu : unchecked 
     242      IF( neuler /= 0 .OR. kt /= nit000 ) THEN           ! remove the forcing from the Asselin filter 
     243         zac = atfp * rdttra(1) 
     244         tb(:,:,1) = tb(:,:,1) - zac * ( qns (:,:) - qns_b (:,:) )         ! non solar surface flux 
     245         sb(:,:,1) = sn(:,:,1) - zac * ( emps(:,:) - emps_b(:,:) )         ! surface salt flux 
     246         ! 
     247         IF( ln_traqsr )                                                   ! solar penetrating flux 
     248            DO jk = 1, nksr 
     249                  DO jj = 1, jpj 
     250                     DO ji = 1, jpi 
     251                        IF( ln_sco ) THEN 
     252                           z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) ) 
     253                        ELSEIF( ln_zps .OR. ln_zco ) THEN 
     254                           z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) *   & 
     255                                &   ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) ) 
     256                        ENDIF 
     257                        zt = zt - zfact1 * z_cor_qsr 
     258                     END DO 
     259                  END DO 
    241260      ENDIF 
    242261      ! 
    … …  
    372391         ENDIF 
    373392      ENDIF 
     393!!gm from Matthieu : unchecked 
     394      ! 
     395      IF( neuler /= 0 .OR. kt /= nit000 ) THEN           ! remove the forcing from the Asselin filter 
     396               IF( lk_vvl ) THEN                                        ! Varying levels 
     397                  DO jj = 1, jpj 
     398                     DO ji = 1, jpi 
     399                        ! - ML - modification for varaiance diagnostics 
     400                        zssh1 = tmask(ji,jj,jk) / ( fse3t(ji,jj,jk) + atfp * ( zfse3ta(ji,jj,jk) - 2*fse3t(ji,jj,jk) & 
     401                           &                                                 + fse3tb(ji,jj,jk) ) ) 
     402                        zt = zssh1 * ( tn(ji,jj,jk) +  atfp * ( ta(ji,jj,jk) - 2 * tn(ji,jj,jk) + tb(ji,jj,jk) ) ) 
     403                        zs = zssh1 * ( sn(ji,jj,jk) +  atfp * ( sa(ji,jj,jk) - 2 * sn(ji,jj,jk) + sb(ji,jj,jk) ) ) 
     404                        ! filter correction for global conservation 
     405                        !------------------------------------------ 
     406                        zfact1 = atfp * rdttra(1) * zssh1 
     407                        IF (jk == 1) THEN                        ! remove sbc trend from time filter 
     408                           zt = zt - zfact1 * ( qn(ji,jj) - qb(ji,jj) ) 
     409!!???                           zs = zs - zfact1 * ( - emp( ji,jj) + empb(ji,jj) ) * zsrau * 35.5e0 
     410                        ENDIF 
     411                        ! remove qsr trend from time filter 
     412                        IF (jk <= nksr) THEN 
     413                           IF( ln_sco ) THEN 
     414                              z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) ) 
     415                           ELSEIF( ln_zps .OR. ln_zco ) THEN 
     416                              z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) *   & 
     417                                   &   ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) ) 
     418                           ENDIF 
     419                           zt = zt - zfact1 * z_cor_qsr 
     420                           IF (jk == nksr) qsrb(ji,jj) = qsr(ji,jj) 
     421                        ENDIF 
     422                        ! - ML - end of modification 
     423                        zssh1 = tmask(ji,jj,1) / zfse3ta(ji,jj,jk) 
     424                        tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1 
     425                        sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1 
     426                     END DO 
     427                  END DO 
     428      ENDIF 
     429!!gm end 
    374430      ! 
    375431   END SUBROUTINE tra_nxt_vvl 

  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90

    r1756 r1870  
    1010   !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate 
    1111   !!            3.2  !  2009-04  (G. Madec & NEMO team)  
     12   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
    1213   !!---------------------------------------------------------------------- 
    1314 
    … …  
    4445   REAL(wp), PUBLIC ::   rn_si2     = 61.8_wp   !: deepest depth of extinction (blue & 0.01 mg.m-3)     (RGB) 
    4546    
    46    ! Module variables 
    4747   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read) 
    48    INTEGER ::   nksr   ! levels below which the light cannot penetrate ( depth larger than 391 m) 
    49    REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption 
     48    
     49   INTEGER                   ::   nksr    ! levels below which the light cannot penetrate (depth larger than 391 m) 
     50   REAL(wp), DIMENSION(3,61) ::   rkrgb   ! tabulated attenuation coefficients for RGB absorption 
    5051 
    5152   !! * Substitutions 
    … …  
    5354#  include "vectopt_loop_substitute.h90" 
    5455   !!---------------------------------------------------------------------- 
    55    !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     56   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    5657   !! $Id$ 
    5758   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    9697      REAL(wp) ::   zchl, zcoef, zsi0r   ! temporary scalars 
    9798      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         - 
     99      REAL(wp) ::   zqsr                 !    -         - 
    98100      REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr            ! 2D workspace 
    99101      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0, ze1 , ze2, ze3, zea    ! 3D workspace 
    … …  
    105107         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    106108         CALL tra_qsr_init 
    107          IF( .NOT.ln_traqsr )   RETURN 
     109         IF( .NOT.ln_traqsr )   RETURN             !!gm not authirized by Coding rules 
    108110      ENDIF 
    109111 
    … …  
    120122            DO jj = 2, jpjm1 
    121123               DO ji = fs_2, fs_jpim1   ! vector opt. 
     124!!gm  how to stecify the mean of time step here : TOP versus OPA time stepping strategy not obvious 
    122125                  ta(ji,jj,jk) = ta(ji,jj,jk) + ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk)  
    123126               END DO 
    … …  
    135138            IF( nn_chldta ==1 ) THEN                             !*  Variable Chlorophyll 
    136139               ! 
    137                CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step 
     140               CALL fld_read( kt, 1, sf_chl )                        ! Read Chl data and provides it at the current time step 
    138141               !          
    139142!CDIR COLLAPSE 
    140143!CDIR NOVERRCHK 
    141                DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl 
     144               DO jj = 1, jpj                                        ! Separation in R-G-B depending of the surface Chl 
    142145!CDIR NOVERRCHK 
    143146                  DO ji = 1, jpi 
    … …  
    149152                  END DO 
    150153               END DO 
    151                ! 
     154               !                                                     ! surface value 
    152155               zsi0r = 1.e0 / rn_si0 
    153                zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B 
    154                ze0(:,:,1) = rn_abs  * qsr(:,:) 
    155                ze1(:,:,1) = zcoef * qsr(:,:) 
    156                ze2(:,:,1) = zcoef * qsr(:,:) 
    157                ze3(:,:,1) = zcoef * qsr(:,:) 
    158                zea(:,:,1) =         qsr(:,:) 
    159                ! 
    160                DO jk = 2, nksr+1 
     156               zcoef = ( 1. - rn_abs ) / 3.e0                           ! equi-partition in R-G-B 
     157!!gm bug !!! zqsr only a constant not an array 
     158               zqsr  = 0.5 * ( qsr_b(:,:) + qsr(:,:) )                  ! mean over 2 time steps 
     159               ze0(:,:,1) = rn_abs  * zqsr 
     160               ze1(:,:,1) = zcoef   * zqsr 
     161               ze2(:,:,1) = zcoef   * zqsr 
     162               ze3(:,:,1) = zcoef   * zqsr 
     163               zea(:,:,1) =           zqsr 
     164               ! 
     165               DO jk = 2, nksr+1                                     ! deeper values 
    161166!CDIR NOVERRCHK 
    162167                  DO jj = 1, jpj 
    163168!CDIR NOVERRCHK    
    164169                     DO ji = 1, jpi 
    165                         zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zsi0r     ) 
     170                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zsi0r       ) 
    166171                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) ) 
    167172                        zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) ) 
    … …  
    175180                  END DO 
    176181               END DO 
    177                ! 
    178                DO jk = 1, nksr                                        ! compute and add qsr trend to ta 
     182!!gm add here  etot3(:,:,jk) = ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk) 
     183               ! 
     184               DO jk = 1, nksr                                       ! compute and add qsr trend to ta 
    179185                  ta(:,:,jk) = ta(:,:,jk) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk) 
    180186               END DO 
    … …  
    184190            ELSE                                                 !*  Constant Chlorophyll 
    185191               DO jk = 1, nksr 
    186                   ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * qsr(:,:) 
     192                  ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) ) 
    187193               END DO 
    188194            ENDIF 
    189  
     195            ! 
    190196         ENDIF 
    191197         !                                                ! ------------------------- ! 
    192198         IF( ln_qsr_2bd ) THEN                            !  2 band light penetration ! 
    193199            !                                             ! ------------------------- ! 
    194             ! 
    195200            DO jk = 1, nksr 
    196201               DO jj = 2, jpjm1 
    197202                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    198                      ta(ji,jj,jk) = ta(ji,jj,jk) + etot3(ji,jj,jk) * qsr(ji,jj) 
     203                     ta(ji,jj,jk) = ta(ji,jj,jk) + etot3(ji,jj,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) ) 
    199204                  END DO 
    200205               END DO 

  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90

    r1739 r1870  
    44   !! Ocean active tracers:  surface boundary condition 
    55   !!============================================================================== 
    6    !! History :  8.2  !  98-10  (G. Madec, G. Roullet, M. Imbard)  Original code 
    7    !!            8.2  !  01-02  (D. Ludicone)  sea ice and free surface 
    8    !!            8.5  !  02-06  (G. Madec)  F90: Free form and module 
     6   !! History :  OPA  !  2998-10  (G. Madec, G. Roullet, M. Imbard)  Original code 
     7   !!            8.2  !  2001-02  (D. Ludicone)  sea ice and free surface 
     8   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module 
     9   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
    910   !!---------------------------------------------------------------------- 
    1011 
    … …  
    3132#  include "vectopt_loop_substitute.h90" 
    3233   !!---------------------------------------------------------------------- 
    33    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     34   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
    3435   !! $Id$ 
    3536   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    4647      !!      and add it to the general trend of tracer equations. 
    4748      !! 
    48       !! ** Method : 
    49       !!      Following Roullet and Madec (2000), the air-sea flux can be divided  
    50       !!      into three effects: (1) Fext, external forcing;  
     49      !! ** Method  :   Following Roullet and Madec (2000), the air-sea flux  
     50      !!      can be divided into three effects:  
     51      !!      (1) Fext, external forcing;  
    5152      !!      (2) Fwi, concentration/dilution effect due to water exchanged  
    5253      !!         at the surface by evaporation, precipitations and runoff (E-P-R);  
    … …  
    5859      !!         solar part is added in traqsr.F routine. 
    5960      !!            ta = ta + q /(rau0 rcp e3t)  for k=1 
    60       !!            - salinity    : no salt flux 
     61      !!            - salinity    : only salt flux exchanged with sea-ice 
    6162      !! 
    6263      !!      The formulation for Fwb and Fwi vary according to the free  
    … …  
    123124      ENDIF 
    124125 
    125       IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration 
     126!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration 
     127      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration 
     128         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns 
     129         qsr(:,:) = 0.e0                     ! qsr set to zero 
     130         IF( kt == nit000 ) THEN             ! idem on before field at nit000 
     131            qns_b(:,:) = qns_b(:,:) + qsr_b(:,:)   
     132            qsr_b(:,:) = 0.e0                      
     133         ENDIF 
     134      ENDIF 
    126135 
    127136      ! Concentration dillution effect on (t,s) 
    128       DO jj = 2, jpj 
    129          DO ji = fs_2, fs_jpim1   ! vector opt. 
    130 #if ! defined key_zco 
    131             zse3t = 1. / fse3t(ji,jj,1) 
    132 #endif 
    133             IF( lk_vvl) THEN 
    134                zta = ro0cpr * qns(ji,jj) * zse3t &                   ! temperature : heat flux 
    135                 &    - emp(ji,jj) * zsrau * tn(ji,jj,1)  * zse3t     ! & cooling/heating effet of EMP flux 
    136                zsa = 0.e0                                            ! No salinity concent./dilut. effect 
    137             ELSE 
    138                zta = ro0cpr * qns(ji,jj) * zse3t     ! temperature : heat flux 
    139                zsa = emps(ji,jj) * zsrau * sn(ji,jj,1)   * zse3t     ! salinity :  concent./dilut. effect 
    140             ENDIF 
    141             ta(ji,jj,1) = ta(ji,jj,1) + zta                          ! add the trend to the general tracer trend 
    142             sa(ji,jj,1) = sa(ji,jj,1) + zsa 
     137!!      DO jj = 2, jpj 
     138!!         DO ji = fs_2, fs_jpim1   ! vector opt. 
     139!!#if ! defined key_zco 
     140!!            zse3t = 1. / fse3t(ji,jj,1) 
     141!!#endif 
     142!!            IF( lk_vvl) THEN 
     143!!               zta = ro0cpr * qns(ji,jj) * zse3t &                   ! temperature : heat flux 
     144!!                &    - emp(ji,jj) * zsrau * tn(ji,jj,1)  * zse3t     ! & cooling/heating effet of EMP flux 
     145!!               zsa = 0.e0                                            ! No salinity concent./dilut. effect 
     146!!            ELSE 
     147!!               zta = ro0cpr * qns(ji,jj) * zse3t     ! temperature : heat flux 
     148!!               zsa = emps(ji,jj) * zsrau * sn(ji,jj,1)   * zse3t     ! salinity :  concent./dilut. effect 
     149!!            ENDIF 
     150!!            ta(ji,jj,1) = ta(ji,jj,1) + zta                          ! add the trend to the general tracer trend 
     151!!            sa(ji,jj,1) = sa(ji,jj,1) + zsa 
     152!!         END DO 
     153!!      END DO 
     154 
     155      !                             ! ---------------------- ! 
     156      IF( lk_vvl ) THEN             !  Variable Volume case  ! 
     157         !                          ! ---------------------- ! 
     158         DO jj = 2, jpj 
     159            DO ji = fs_2, fs_jpim1   ! vector opt. 
     160               zse3t = 0.5 / fse3t(ji,jj,1) 
     161               !                           ! temperature: heat flux + cooling/heating effet of EMP flux 
     162               ta(ji,jj,1) = ta(ji,jj,1) + (  ro0cpr * ( qns(ji,jj) + qns_b(ji,jj) )                &  
     163                  &                         - zsrau  * ( emp(ji,jj) + emp_b(ji,jj) ) * tn(ji,jj,1)  ) * zse3t 
     164               !                           ! salinity: salt flux 
     165               sa(ji,jj,1) = sa(ji,jj,1) + ( emps(ji,jj) + emps_b(ji,jj) ) * zse3t 
     166 
     167!!gm BUG : in key_vvl emps must be modified to only include the salt flux due to with sea-ice freezing/melting 
     168!!gm       otherwise this flux will be missing  ==> modification required in limsbc,  limsbc_2 and CICE interface. 
     169 
     170            END DO 
    143171         END DO 
    144       END DO 
     172         !                          ! ---------------------- ! 
     173      ELSE                          !  Constant Volume case  ! 
     174         !                          ! ---------------------- ! 
     175         DO jj = 2, jpj 
     176            DO ji = fs_2, fs_jpim1   ! vector opt. 
     177               zse3t = 0.5 / fse3t(ji,jj,1) 
     178               !                           ! temperature: heat flux  
     179               ta(ji,jj,1) = ta(ji,jj,1) + ro0cpr * ( qns (ji,jj) + qns_b (ji,jj) )               * zse3t  
     180               !                           ! salinity: salt flux + concent./dilut. effect (both in emps) 
     181               sa(ji,jj,1) = sa(ji,jj,1) + zsrau  * ( emps(ji,jj) + emps_b(ji,jj) ) * sn(ji,jj,1) * zse3t 
     182            END DO 
     183         END DO 
     184         ! 
     185      ENDIF 
    145186 
    146       IF( l_trdtra ) THEN      ! save the sbc trends for diagnostic 
     187      IF( l_trdtra ) THEN        ! save the sbc trends for diagnostic 
    147188         ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:) 
    148189         ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:) 
    149190         CALL trd_mod(ztrdt, ztrds, jptra_trd_nsr, 'TRA', kt) 
    150191      ENDIF 
    151       ! 
     192      !                          ! control print 
    152193      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' sbc  - Ta: ', mask1=tmask,   & 
    153194         &                       tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
Note: See TracChangeset for help on using the changeset viewer.