Changeset 1870


Ignore:
Timestamp:
2010-05-12T17:36:00+02:00 (10 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.