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

Changeset 700


Ignore:
Timestamp:
2007-10-09T14:24:27+02:00 (17 years ago)
Author:
smasson
Message:

sbc(1): bugfix on wind stress computation #2

Location:
trunk/NEMO
Files:
2 added
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/LIM_SRC/icestp.F90

    r699 r700  
    2121   USE flx_oce         ! forcings variables 
    2222   USE dom_ice 
    23    USE cpl_oce 
    2423   USE daymod 
    2524   USE phycst          ! Define parameters for the routines 
     
    4847   !!----------------------------------------------------- 
    4948   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    50    !! $Id$ 
     49   !! $Header$  
    5150   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5251   !!----------------------------------------------------- 
     
    7170 
    7271      INTEGER                      ::   ji, jj   ! dummy loop indices 
    73       REAL(wp), DIMENSION(jpi,jpj) ::   zsss_io, zsss2_io, zsss3_io          ! tempory workspaces 
     72      REAL(wp), DIMENSION(jpi,jpj) ::   zsss_io, zsss2_io, zsss3_io   ! tempory workspaces 
     73      REAL(wp)                     ::   zmult                         ! tempory workspaces 
    7474      !!---------------------------------------------------------------------- 
    7575 
    7676      IF( kt == nit000 ) THEN 
    77          IF( lk_cpl ) THEN 
    78             IF(lwp) WRITE(numout,*) 
    79             IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)' 
    80             IF(lwp) WRITE(numout,*) '~~~~~~~   coupled case' 
    81          ELSE 
    82             IF(lwp) WRITE(numout,*) 
    83             IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)'  
    84             IF(lwp) WRITE(numout,*) '~~~~~~~   forced case using bulk formulea' 
    85          ENDIF 
    86          !  Initialize fluxes fields 
     77         IF(lwp) WRITE(numout,*) 
     78         IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)' 
     79          !  Initialize fluxes fields 
    8780         gtaux(:,:) = 0.e0 
    8881         gtauy(:,:) = 0.e0 
     
    9790      !    cumulate fields 
    9891      ! 
    99       sst_io(:,:) = sst_io(:,:) + tn(:,:,1) + rt0 
     92      sst_io(:,:) = sst_io(:,:) + tn(:,:,1) 
    10093      sss_io(:,:) = sss_io(:,:) + sn(:,:,1) 
    10194 
    10295  
    103       ! vectors at F-point 
     96      ! surface ocean current at I-point (from U-, V-points) 
    10497      DO jj = 2, jpj 
    10598         DO ji = fs_2, jpi   ! vector opt. 
    106             u_io(ji,jj) = u_io(ji,jj) + 0.5 * ( un(ji-1,jj  ,1) + un(ji-1,jj-1,1) ) 
    107             v_io(ji,jj) = v_io(ji,jj) + 0.5 * ( vn(ji  ,jj-1,1) + vn(ji-1,jj-1,1) ) 
     99            u_io (ji,jj) = u_io(ji,jj) + un(ji-1,jj  ,1) + un(ji-1,jj-1,1) 
     100            v_io (ji,jj) = v_io(ji,jj) + vn(ji  ,jj-1,1) + vn(ji-1,jj-1,1) 
    108101         END DO 
    109102      END DO 
     
    113106          
    114107         ! The LIM model is going to be call 
    115          sst_io(:,:) = sst_io(:,:) / FLOAT( nfice ) * tmask(:,:,1) 
    116          sss_io(:,:) = sss_io(:,:) / FLOAT( nfice ) 
     108         zmult = 1. / REAL( nfice, wp ) 
     109         sst_io(:,:) = zmult * sst_io(:,:) * tmask(:,:,1) + rt0 
     110         sss_io(:,:) = zmult * sss_io(:,:) 
    117111 
    118112         ! stress from ocean U- and V-points to ice U,V point 
     113         zmult = 0.5 * zmult 
    119114         DO jj = 2, jpj 
    120115            DO ji = fs_2, jpi   ! vector opt. 
    121                gtaux(ji,jj) = 0.5 * ( taux(ji-1,jj  ) + taux(ji-1,jj-1) ) 
    122                gtauy(ji,jj) = 0.5 * ( tauy(ji  ,jj-1) + tauy(ji-1,jj-1) ) 
    123                u_io  (ji,jj) = u_io(ji,jj) / FLOAT( nfice ) 
    124                v_io  (ji,jj) = v_io(ji,jj) / FLOAT( nfice ) 
     116               ! ice/atmosphere stress at I-point (from U-, V-points) 
     117               gtaux(ji,jj) = 0.5 * ( taux_ice(ji-1,jj  ) + taux_ice(ji-1,jj-1) )  
     118               gtauy(ji,jj) = 0.5 * ( tauy_ice(ji  ,jj-1) + tauy_ice(ji-1,jj-1) ) 
     119                
     120               u_io (ji,jj) = zmult * u_io(ji,jj) 
     121               v_io (ji,jj) = zmult * v_io(ji,jj) 
    125122            END DO 
    126123         END DO 
     
    168165!!gm      END DO 
    169166!!gm 
    170       zsss_io (:,:) = SQRT( sss_io(:,:) )  
    171       zsss2_io(:,:) =  sss_io(:,:) *  sss_io(:,:) 
    172       zsss3_io(:,:) = zsss_io(:,:) * zsss_io(:,:) * zsss_io(:,:) 
    173  
    174       DO jj = 1, jpj 
    175          DO ji = 1, jpi 
    176             tfu(ji,jj)  = ABS ( rt0 - 0.0575       *   sss_io(ji,jj)   & 
    177                &                    + 1.710523e-03 * zsss3_io(ji,jj)   & 
    178                &                    - 2.154996e-04 * zsss2_io(ji,jj) ) * tms(ji,jj) 
     167         zsss_io (:,:) = SQRT( sss_io(:,:) )  
     168         zsss2_io(:,:) =  sss_io(:,:) *  sss_io(:,:) 
     169         zsss3_io(:,:) = zsss_io(:,:) * zsss_io(:,:) * zsss_io(:,:) 
     170          
     171         DO jj = 1, jpj 
     172            DO ji = 1, jpi 
     173               tfu(ji,jj)  = ABS ( rt0 - 0.0575       *   sss_io(ji,jj)   & 
     174                    &                    + 1.710523e-03 * zsss3_io(ji,jj)   & 
     175                    &                    - 2.154996e-04 * zsss2_io(ji,jj) ) * tms(ji,jj) 
     176            END DO 
    179177         END DO 
    180       END DO 
    181     
    182178          
    183179 
     
    196192 
    197193         !                                                           !-----------------------! 
    198          CALL lim_rst_opn( kt )                                   ! Open Ice restart file ! 
     194         CALL lim_rst_opn( kt )                                      ! Open Ice restart file ! 
    199195         !                                                           !-----------------------! 
    200196 
     
    240236 
    241237         IF( MOD( kt + nfice -1, ninfo ) == 0 .OR. ntmoy == 1 )  THEN        !-----------------! 
    242             CALL lim_dia( kt )                                    ! Ice Diagnostics ! 
     238            CALL lim_dia( kt )                                       ! Ice Diagnostics ! 
    243239         ENDIF                                                       !-----------------! 
    244240 
     
    248244 
    249245         !                                                           !------------------------! 
    250          IF( lrst_ice ) CALL lim_rst_write( kt )                  ! Write Ice restart file ! 
     246         IF( lrst_ice ) CALL lim_rst_write( kt )                     ! Write Ice restart file ! 
    251247         !                                                           !------------------------! 
    252248 
     
    262258         fr2_i0  (:,:) = 0.e0 
    263259         evap    (:,:) = 0.e0 
    264 #if defined key_coupled  
    265          rrunoff (:,:) = 0.e0 
    266          calving (:,:) = 0.e0 
    267 #else 
     260#if ! defined key_coupled  
    268261         qla_ice (:,:) = 0.e0 
    269262         dqla_ice(:,:) = 0.e0 
  • trunk/NEMO/LIM_SRC/limflx.F90

    r699 r700  
    3939   !!---------------------------------------------------------------------- 
    4040   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    41    !! $Id$ 
     41   !! $Header$  
    4242   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    4343   !!---------------------------------------------------------------------- 
     
    106106      DO jj = 1, jpj 
    107107         DO ji = 1, jpi 
    108             zinda   = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 
    109             ifvt    = zinda  *  MAX( rzero , SIGN( rone, -phicif(ji,jj) ) ) 
    110             i1mfr   = 1.0 - MAX( rzero , SIGN( rone ,  - ( 1.0 - frld(ji,jj) ) ) ) 
    111             idfr    = 1.0 - MAX( rzero , SIGN( rone , frld(ji,jj) - pfrld(ji,jj) ) ) 
    112             iflt    = zinda  * (1 - i1mfr) * (1 - ifvt ) 
    113             ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr 
    114             iadv    = ( 1  - i1mfr ) * zinda 
    115             ifral   = ( 1  - i1mfr * ( 1 - ial ) )    
     108!!sm            zinda   = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 
     109!!sm            ifvt    = zinda  *  MAX( rzero , SIGN( rone, -phicif(ji,jj) ) ) 
     110!!sm            i1mfr   = 1.0 - MAX( rzero , SIGN( rone ,  - ( 1.0 - frld(ji,jj) ) ) ) 
     111!!sm            idfr    = 1.0 - MAX( rzero , SIGN( rone , frld(ji,jj) - pfrld(ji,jj) ) ) 
     112 
     113            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   !   = 0. if pure ocean else 1. (at previous time) 
     114 
     115            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   !   = 0. if pure ocean else 1. (at current  time) 
     116 
     117            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      !   = 1. if (snow and no ice at previous time) else 0. ??? 
     118            ELSE                             ;   ifvt = 0. 
     119            ENDIF 
     120 
     121            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases from previous to current 
     122            ELSE                                     ;   idfr = 1.    
     123            ENDIF 
     124 
     125            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous and pure ocean at current 
     126 
     127            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr 
     128!                 snow no ice   ice         ice or nothing  lead fraction increases 
     129!                 at previous   now           at previous 
     130!                -> ice aera increases  ???         -> ice aera decreases ??? 
     131 
     132            iadv    = ( 1  - i1mfr ) * zinda   
     133!                     pure ocean      ice at 
     134!                     at current      previous 
     135!                        -> = 1. if ice disapear between previous and current 
     136 
     137            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
     138!                            ice at     ??? 
     139!                            current          
     140!                         -> ??? 
     141  
    116142            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv  
     143!                                                    ice disapear                            
     144! 
     145! 
    117146            z1mthcm =   1. - thcm(ji,jj)        
    118147            !   computation the solar flux at ocean surface 
     
    179208      !-----------------------------------------------! 
    180209 
    181       ftaux (:,:) = - tio_u(:,:) * rau0   ! taux ( ice: N/m2/rau0, ocean: N/m2 ) 
    182       ftauy (:,:) = - tio_v(:,:) * rau0   ! tauy ( ice: N/m2/rau0, ocean: N/m2 )                 
    183210      freeze(:,:) = 1.0 - frld(:,:)       ! Sea ice cover             
    184211      tn_ice(:,:) = sist(:,:)             ! Ice surface temperature                       
     
    201228         CALL prt_ctl(tab2d_1=fsolar, clinfo1=' lim_flx: fsolar : ', tab2d_2=fnsolar, clinfo2=' fnsolar : ') 
    202229         CALL prt_ctl(tab2d_1=fmass , clinfo1=' lim_flx: fmass  : ', tab2d_2=fsalt  , clinfo2=' fsalt   : ') 
    203          CALL prt_ctl(tab2d_1=ftaux , clinfo1=' lim_flx: ftaux  : ', tab2d_2=ftauy  , clinfo2=' ftauy   : ') 
    204230         CALL prt_ctl(tab2d_1=freeze, clinfo1=' lim_flx: freeze : ', tab2d_2=tn_ice , clinfo2=' tn_ice  : ') 
    205231      ENDIF  
  • trunk/NEMO/LIM_SRC/limhdf.F90

    r699 r700  
    3434   !!---------------------------------------------------------------------- 
    3535   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    36    !! $Id$ 
     36   !! $Header$  
    3737   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    3838   !!---------------------------------------------------------------------- 
  • trunk/NEMO/LIM_SRC/limrhg.F90

    r699 r700  
    44   !!   Ice rheology :  performs sea ice rheology 
    55   !!====================================================================== 
     6   !! History :   0.0  !  93-12  (M.A. Morales Maqueda.)  Original code 
     7   !!             1.0  !  94-12  (H. Goosse)  
     8   !!             2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp 
     9   !!             3.0  !  07-06  (S. Masson, G. Madec)  ice/atmosphere stress 
     10   !!                             provided at I-point in forced and coupled mode 
     11   !!---------------------------------------------------------------------- 
    612#if defined key_ice_lim 
    713   !!---------------------------------------------------------------------- 
     
    3339   !!---------------------------------------------------------------------- 
    3440   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    35    !! $Id$ 
     41   !! $Header$  
    3642   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    3743   !!---------------------------------------------------------------------- 
     
    5056      !! ** Action  : - compute u_ice, v_ice the sea-ice velocity 
    5157      !! 
    52       !! History : 
    53       !!   0.0  !  93-12  (M.A. Morales Maqueda.)  Original code 
    54       !!   1.0  !  94-12  (H. Goosse)  
    55       !!   2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp 
    5658      !!------------------------------------------------------------------- 
    57       ! * Arguments 
    5859      INTEGER, INTENT(in) ::   k_j1 ,  &  ! southern j-index for ice computation 
    5960         &                     k_jpj      ! northern j-index for ice computation 
    6061 
    61       ! * Local variables 
    6262      INTEGER ::   ji, jj              ! dummy loop indices 
    6363 
     
    7777         zresm, zunw, zvnw, zur, zvr, zmod, za, zac, & 
    7878         zmpzas, zstms, zindu, zindu1, zusdtp, zmassdt, zcorlal,  & 
    79          ztrace2, zdeter, zdelta, zsang, zmask, zdgp, zdgi, zdiag 
     79         ztrace2, zdeter, zdelta, zsang, zmask, zdgp, zdgi, zdiag, zic 
    8080      REAL(wp),DIMENSION(jpi,jpj) ::   & 
    8181         zpresh, zfrld, zmass, zcorl,     & 
     
    9292      zv0(:,:) = v_ice(:,:) 
    9393 
    94       ! Ice mass, ice strength, and wind stress at the center            | 
    95       ! of the grid squares.                                             | 
     94      ! masked Ice mass, ice strength, Ice-cover fraction and Lead faction  at T-point  
    9695      !------------------------------------------------------------------- 
    9796 
     
    10099            za1(ji,jj)    = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) ) 
    101100            zpresh(ji,jj) = tms(ji,jj) *  pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) ) 
    102 #if defined key_lim_cp1 && defined key_coupled 
    103             zb1(ji,jj)    = tms(ji,jj) * gtaux(ji,jj) * ( 1.0 - frld(ji,jj) ) 
    104             zb2(ji,jj)    = tms(ji,jj) * gtauy(ji,jj) * ( 1.0 - frld(ji,jj) ) 
    105 #else 
    106             zb1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 
    107             zb2(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 
    108 #endif 
     101 
     102            zb1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )     
     103            zb2(ji,jj)    = tms(ji,jj) *         frld(ji,jj)       
    109104         END DO 
    110105      END DO 
    111106 
    112  
    113       !--------------------------------------------------------------------------- 
    114       !  Wind stress, coriolis and mass terms at the corners of the grid squares | 
    115       !  Gradient of ice strenght.                                               | 
     107      !  Wind stress, coriolis, mass and Gradient of ice strenght at I-point 
    116108      !--------------------------------------------------------------------------- 
    117109          
     
    122114            zusw  = 1.0 / MAX( zstms, epsd ) 
    123115 
    124             zt11 = tms(ji  ,jj  ) * frld(ji  ,jj  )  
    125             zt12 = tms(ji-1,jj  ) * frld(ji-1,jj  )  
    126             zt21 = tms(ji  ,jj-1) * frld(ji  ,jj-1)  
    127             zt22 = tms(ji-1,jj-1) * frld(ji-1,jj-1) 
    128  
    129             ! Leads area. 
    130             zfrld(ji,jj) =  (  zt11 * wght(ji,jj,2,2) + zt12 * wght(ji,jj,1,2)   & 
    131                &             + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw 
     116            ! Leads area at I-point 
     117            zfrld(ji,jj) = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     118               &           + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
     119 
     120            ! Ice cover area at I-point 
     121            zic          = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     122               &           + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
     123 
     124            ! Wind stress. 
     125            ztagnx = zic * gtaux(ji,jj) 
     126            ztagny = zic * gtauy(ji,jj) 
    132127 
    133128            ! Mass and coriolis coeff. 
     
    135130               &           + za1(ji,jj-1) * wght(ji,jj,2,1) + za1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
    136131            zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj) 
    137  
    138             ! Wind stress. 
    139 #if defined key_lim_cp1 && defined key_coupled 
    140             ztagnx = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    141                &     + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
    142             ztagny = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    143                &     + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
    144 #else 
    145             ztagnx = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    146                &     + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtaux(ji,jj) 
    147             ztagny = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    148                &     + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtauy(ji,jj) 
    149 #endif 
    150132 
    151133            ! Gradient of ice strength 
     
    162144            ! Computation of the velocity field taking into account the ice-ice interaction.                                  
    163145            ! Terms that are independent of the velocity field. 
    164             za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v_oce(ji,jj) - zgphsx 
    165             za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u_oce(ji,jj) - zgphsy 
     146            za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v_io(ji,jj) - zgphsx 
     147            za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u_io(ji,jj) - zgphsy 
    166148         END DO 
    167149      END DO 
     
    403385               zsang  = SIGN( 1.e0, fcor(1,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere 
    404386               DO ji = 2, jpim1 
    405                  zur     = u_ice(ji,jj) - u_oce(ji,jj) 
    406                  zvr     = v_ice(ji,jj) - v_oce(ji,jj) 
     387                 zur     = u_ice(ji,jj) - u_io(ji,jj) 
     388                 zvr     = v_ice(ji,jj) - v_io(ji,jj) 
    407389                 zmod    = SQRT( zur * zur + zvr * zvr) * ( 1.0 - zfrld(ji,jj) ) 
    408390                 za      = rhoco * zmod 
     
    413395 
    414396                  za1(ji,jj) =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   & 
    415                      &        + za * ( cangvg * u_oce(ji,jj) - zsang * v_oce(ji,jj) ) 
     397                     &        + za * ( cangvg * u_io(ji,jj) - zsang * v_io(ji,jj) ) 
    416398 
    417399                  za2(ji,jj) =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   & 
    418                      &        + za * ( cangvg * v_oce(ji,jj) + zsang * u_oce(ji,jj) ) 
     400                     &        + za * ( cangvg * v_io(ji,jj) + zsang * u_io(ji,jj) ) 
    419401 
    420402                  zb1(ji,jj)  = zmassdt + zac - zc1u(ji,jj) 
  • trunk/NEMO/OPA_SRC/SBC/ocesbc.F90

    r699 r700  
    1212   USE oce            ! dynamics and tracers variables 
    1313   USE dom_oce        ! ocean space domain variables 
    14    USE cpl_oce        ! coupled ocean-atmosphere variables 
    1514   USE ice_oce        ! sea-ice variable 
    1615   USE blk_oce        ! bulk variables 
     
    2827   USE in_out_manager ! I/O manager 
    2928   USE prtctl         ! Print control 
     29   USE diurnalcycle   !  
    3030 
    3131   IMPLICIT NONE 
     
    3838   REAL(wp), PUBLIC ::   &  !: 
    3939      aplus, aminus,     &  !: 
    40       empold = 0.e0         !: current year freshwater budget correction 
     40      empold = 0.e0,     &  !: current year freshwater budget correction 
     41      area_g                !: surface og the global ocean 
    4142   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 
    42       qt  ,         &  !: total surface heat flux (w/m2) 
    43       qsr ,         &  !: solar radiation (w/m2) 
     43!CT      qt  ,         &  !: total surface heat flux (w/m2) 
     44!CT      qsr ,         &  !: solar radiation (w/m2) 
    4445      qla ,         &  !: Latent heat flux (W/m2) 
    4546      qlw ,         &  !: Long Wave Heat Flux (W/m2) 
     
    5253   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 
    5354      dmp              !: internal dampind term 
     55   REAL(wp), DIMENSION(jpi,jpj) ::   &  !: 
     56      e1e2_i    ! area of the interior domain (e1t*e2t*tmask_i) 
    5457#endif 
    5558 
     
    5962   !!---------------------------------------------------------------------- 
    6063   !!  OPA 9.0 , LOCEAN-IPSL (2005)  
    61    !! $Id$ 
     64   !! $Header$  
    6265   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    6366   !!---------------------------------------------------------------------- 
     
    9093      !! * Local declarations 
    9194      INTEGER ::   ji, jj                   ! dummy loop indices 
    92       REAL(wp) ::   ztx, ztaux, zty, ztauy 
    93       REAL(wp) ::   ztdta, ztgel, zqrp 
    94       !!---------------------------------------------------------------------- 
     95      REAL(wp) ::   ztx, ztaux_ice, zty, ztauy_ice 
     96      REAL(wp) ::   ztdta, ztgel, zqrp, zm_emp 
     97      !!---------------------------------------------------------------------- 
     98 
     99      IF( kt == nit000 ) THEN 
     100         e1e2_i(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:) 
     101         area_g = SUM( e1e2_i(:,:) ) 
     102         IF( lk_mpp )   CALL  mpp_sum( area_g )   ! sum over the global domain 
     103      ENDIF 
    95104 
    96105      ! 1. initialization to zero at kt = nit000 
    97106      ! --------------------------------------- 
    98107 
    99       IF( kt == nit000 ) THEN      
    100          qsr   (:,:) = 0.e0 
    101          freeze(:,:) = 0.e0 
    102          qt    (:,:) = 0.e0 
    103          qrp   (:,:) = 0.e0 
    104          emp   (:,:) = 0.e0 
    105          emps  (:,:) = 0.e0 
    106          erp   (:,:) = 0.e0 
    107 #if ! defined key_dynspg_rl  
    108          dmp   (:,:) = 0.e0 
    109 #endif 
    110       ENDIF 
    111  
    112108      IF( MOD( kt-1, nfice ) == 0 ) THEN  
    113109 
    114          CALL oce_sbc_dmp   ! Computation of internal and evaporation damping terms        
    115  
    116          ! Surface heat flux (W/m2) 
    117          ! ----------------------- 
    118  
    119          ! restoring heat flux 
    120          DO jj = 1, jpj 
    121             DO ji = 1, jpi 
    122                ztgel = fzptn(ji,jj) 
    123 #if defined key_dtasst 
    124                ztdta = MAX( sst(ji,jj),    ztgel ) 
    125 #else 
    126                ztdta = MAX( t_dta(ji,jj,1), ztgel ) 
    127 #endif 
    128                zqrp = dqdt0 * ( tb(ji,jj,1) - ztdta ) 
    129  
    130                qrp(ji,jj) = (1.0-freeze(ji,jj) ) * zqrp 
    131             END DO 
    132          END DO 
    133  
    134110         ! non solar heat flux + solar flux + restoring 
    135          qt (:,:) = fnsolar(:,:) + fsolar(:,:) + qrp(:,:) 
     111         qt (:,:) = fnsolar(:,:) + fsolar(:,:) 
    136112 
    137113         ! solar flux 
     
    140116#if ! defined key_dynspg_rl   
    141117         ! total concentration/dilution effect (use on SSS) 
    142          emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:) 
    143  
     118         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + calving(:,:) 
    144119         ! total volume flux (use on sea-surface height) 
    145          emp (:,:) = fmass(:,:)  -  dmp(:,:) + runoff(:,:) + erp(:,:) 
     120         emp (:,:) = fmass(:,:)              + runoff(:,:) + calving(:,:) 
     121 
     122         ! Compute the mean emp(:,:) 
     123         zm_emp =  SUM( emp(:,:)*e1e2_i(:,:) ) / area_g 
     124         IF( lk_mpp )   CALL  mpp_sum( zm_emp )   ! sum over the global domain 
     125 
     126         ! Correct emp() & emps() in substracting the emp() mean 
     127         emps(:,:) = emps(:,:) - zm_emp 
     128         emp (:,:) = emp (:,:) - zm_emp 
    146129#else 
    147130         ! Rigid-lid (emp=emps=E-P-R+Erp)  
    148131         ! freshwater flux 
    149          emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:) 
     132         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + calving(:,:) 
    150133         emp (:,:) = emps(:,:) 
    151134#endif 
     
    153136         DO jj = 1, jpjm1 
    154137            DO ji = 1, fs_jpim1   ! vertor opt. 
    155                ztx   = 0.5 * ( freeze(ji+1,jj) + freeze(ji+1,jj+1) ) 
    156                ztaux = 0.5 * ( ftaux (ji+1,jj) + ftaux (ji+1,jj+1) ) 
    157                taux(ji,jj) = (1.0-ztx) * taux(ji,jj) + ztx * ztaux 
    158  
    159                zty   = 0.5 * ( freeze(ji,jj+1) + freeze(ji+1,jj+1) ) 
    160                ztauy = 0.5 * ( ftauy (ji,jj+1) + ftauy (ji+1,jj+1) ) 
    161                tauy(ji,jj) = (1.0-zty) * tauy(ji,jj) + zty * ztauy 
     138               ztx         = 0.5 * ( freeze(ji  ,jj) + freeze(ji+1,jj  ) )   !   ice cover        at U-point 
     139               ztaux_ice   = 0.5 * ( ftaux (ji+1,jj) + ftaux (ji+1,jj+1) )   !   ice/ocean stress at U-point 
     140               taux(ji,jj) = (1.0-ztx) * taux_oce(ji,jj) + ztx * ztaux_ice 
     141 
     142               zty         = 0.5 * ( freeze(ji,jj  ) + freeze(ji  ,jj+1) )   !   ice cover        at V-point 
     143               ztauy_ice   = 0.5 * ( ftauy (ji,jj+1) + ftauy (ji+1,jj+1) )   !   ice/ocean stress at V-point 
     144               tauy(ji,jj) = (1.0-zty) * tauy_oce(ji,jj) + zty * ztauy_ice 
    162145            END DO 
    163146         END DO 
     
    200183      !! * Local declarations 
    201184      INTEGER  ::   ji, jj                   ! dummy loop indices 
    202       REAL(wp) ::   ztx, ztaux, zty, ztauy 
     185      REAL(wp) ::   ztx, ztaux, zty, ztauy, zm_emp 
    203186      !!---------------------------------------------------------------------- 
    204187 
     
    219202         dmp    (:,:) = 0.e0 
    220203#endif 
     204         e1e2_i(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:) 
     205         area_g = SUM( e1e2_i(:,:) ) 
     206 
     207         IF( lk_mpp )   CALL  mpp_sum( area_g )   ! sum over the global domain 
    221208      ENDIF 
    222209#if defined key_flx_core 
     
    244231         ! total volume flux (use on sea-surface height) 
    245232         emp (:,:) = fmass(:,:) -   dmp(:,:) + runoff(:,:) + erp(:,:) + empold       
     233 
     234         ! Compute the mean emp(:,:) 
     235         zm_emp =  SUM( emp(:,:)*e1e2_i(:,:) ) / area_g 
     236 
     237         IF( lk_mpp )   CALL  mpp_sum( zm_emp )   ! sum over the global domain 
     238 
     239         ! Correct emp() & emps() in substracting the emp() mean 
     240         emps(:,:) = emps(:,:) - zm_emp 
     241         emp (:,:) = emp (:,:) - zm_emp 
    246242#else 
    247243         ! Rigid-lid (emp=emps=E-P-R+Erp) 
     
    278274      ENDIF 
    279275 
     276      IF( ln_dcy )   CALL diurnal_cycle( kt ) 
     277 
    280278   END SUBROUTINE oce_sbc 
    281279 
     
    307305      !!   2.0  !  02-12  (G. Madec)  F90: Free form and module 
    308306      !!---------------------------------------------------------------------- 
    309       !! * Modules used 
    310       USE cpl_oce       ! coupled ocean-atmosphere variables 
    311  
    312307      !! * Arguments 
    313308      INTEGER, INTENT( in  ) ::   kt   ! ocean time step index 
     
    804799      REAL(wp), DIMENSION(jpi,jpj)  :: zsss, zfreeze 
    805800      REAL(wp) ::   zerp, zsrp 
    806       CHARACTER (len=71) :: charout 
    807 #if ! defined key_dynspg_rl 
    808       REAL(wp) ::   zwei 
    809       REAL(wp) ::   zerpplus(jpi,jpj), zerpminus(jpi,jpj) 
    810       REAL(wp) ::   zplus, zminus, zadefi 
    811 # if defined key_tradmp 
    812       INTEGER jk 
    813       REAL(wp), DIMENSION(jpi,jpj) ::   zstrdmp 
    814 # endif 
    815 #endif 
     801!CT      CHARACTER (len=71) :: charout 
     802!CT#if ! defined key_dynspg_rl 
     803!CT      REAL(wp) ::   zwei 
     804!CT      REAL(wp) ::   zerpplus(jpi,jpj), zerpminus(jpi,jpj) 
     805!CT      REAL(wp) ::   zplus, zminus, zadefi 
     806!CT# if defined key_tradmp 
     807!CT      INTEGER jk 
     808!CT      REAL(wp), DIMENSION(jpi,jpj) ::   zstrdmp 
     809!CT# endif 
     810!CT#endif 
    816811      !!---------------------------------------------------------------------- 
    817812 
     
    839834          
    840835      ! Internal damping 
    841 # if defined key_tradmp 
    842       ! Vertical mean of dampind trend (computed in tradmp module) 
    843       zstrdmp(:,:) = 0.e0 
    844       DO jk = 1, jpk 
    845          zstrdmp(:,:) = zstrdmp(:,:) + strdmp(:,:,jk) * fse3t(:,:,jk) 
    846       END DO 
    847       ! volume flux associated to internal damping to climatology 
    848       dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + 1.e-20 ) 
    849 # else 
     836!CT# if defined key_tradmp 
     837!CT      ! Vertical mean of dampind trend (computed in tradmp module) 
     838!CT      zstrdmp(:,:) = 0.e0 
     839!CT      DO jk = 1, jpk 
     840!CT         zstrdmp(:,:) = zstrdmp(:,:) + strdmp(:,:,jk) * fse3t(:,:,jk) 
     841!CT      END DO 
     842!CT      ! volume flux associated to internal damping to climatology 
     843!CT      dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + 1.e-20 ) 
     844!CT# else 
    850845      dmp(:,:) = 0.e0            ! No internal damping 
    851 # endif 
     846!CT# endif 
    852847       
    853848      !   evaporation damping term ( Surface restoring ) 
    854       zerpplus (:,:) = 0.e0 
    855       zerpminus(:,:) = 0.e0 
    856       zplus  =  15. / rday 
    857       zminus = -15. / rday 
     849!CT      zerpplus (:,:) = 0.e0 
     850!CT      zerpminus(:,:) = 0.e0 
     851!CT      zplus  =  15. / rday 
     852!CT      zminus = -15. / rday 
    858853       
    859854      DO jj = 1, jpj 
     
    863858               & / ( zsss(ji,jj) + 1.e-20        ) 
    864859             
    865             zerp = MIN( zerp, zplus  ) 
    866             zerp = MAX( zerp, zminus ) 
     860!CT            zerp = MIN( zerp, zplus  ) 
     861!CT            zerp = MAX( zerp, zminus ) 
    867862            erp(ji,jj) = zerp 
    868             zerpplus (ji,jj) = MAX( erp(ji,jj), 0.e0 ) 
    869             zerpminus(ji,jj) = MIN( erp(ji,jj), 0.e0 ) 
     863!CT            zerpplus (ji,jj) = MAX( erp(ji,jj), 0.e0 ) 
     864!CT            zerpminus(ji,jj) = MIN( erp(ji,jj), 0.e0 ) 
    870865         END DO 
    871866      END DO 
    872867 
    873       aplus  = 0.e0 
    874       aminus = 0.e0 
    875  
    876       IF( nbit_cmp == 1) THEN 
    877  
    878          IF(ln_ctl)   THEN 
    879             WRITE(charout,FMT="('oce_sbc_dmp : a+ = ',D23.16, ' a- = ',D23.16)") aplus, aminus 
    880             CALL prt_ctl_info(charout) 
    881          ENDIF 
    882          erp(:,:) = 0.e0 
    883  
    884       ELSE 
    885  
    886          DO jj = 1, jpj 
    887             DO ji = 1, jpi 
    888                zwei   = e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj) 
    889                aplus  = aplus  + zerpplus (ji,jj) * zwei 
    890                aminus = aminus - zerpminus(ji,jj) * zwei 
    891             END DO 
    892          END DO 
    893          IF( lk_mpp )   CALL mpp_sum( aplus  )   ! sums over the global domain 
    894          IF( lk_mpp )   CALL mpp_sum( aminus ) 
    895  
    896          IF(ln_ctl)   THEN 
    897             WRITE(charout,FMT="('oce_sbc_dmp : a+ = ',D23.16, ' a- = ',D23.16)") aplus, aminus 
    898             CALL prt_ctl_info(charout) 
    899          ENDIF 
    900  
    901          zadefi = MIN( aplus, aminus ) 
    902          IF( zadefi == 0.e0 ) THEN  
    903             erp(:,:) = 0.e0 
    904          ELSE 
    905             erp(:,:) = zadefi * ( zerpplus(:,:) / aplus + zerpminus(:,:) / aminus ) 
    906          ENDIF 
    907  
    908       END IF 
     868!CT      aplus  = 0.e0 
     869!CT      aminus = 0.e0 
     870!CT       
     871!CT      IF( nbit_cmp == 1) THEN 
     872!CT          
     873!CT         IF(ln_ctl)   THEN 
     874!CT            WRITE(charout,FMT="('oce_sbc_dmp : a+ = ',D23.16, ' a- = ',D23.16)") aplus, aminus 
     875!CT            CALL prt_ctl_info(charout) 
     876!CT         ENDIF 
     877!CT         erp(:,:) = 0.e0 
     878!CT 
     879!CT      ELSE 
     880!CT 
     881!CT         DO jj = 1, jpj 
     882!CT            DO ji = 1, jpi 
     883!CT               zwei   = e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj) 
     884!CT               aplus  = aplus  + zerpplus (ji,jj) * zwei 
     885!CT               aminus = aminus - zerpminus(ji,jj) * zwei 
     886!CT            END DO 
     887!CT         END DO 
     888!CT         IF( lk_mpp )   CALL mpp_sum( aplus  )   ! sums over the global domain 
     889!CT         IF( lk_mpp )   CALL mpp_sum( aminus ) 
     890!CT 
     891!CT         IF(ln_ctl)   THEN 
     892!CT            WRITE(charout,FMT="('oce_sbc_dmp : a+ = ',D23.16, ' a- = ',D23.16)") aplus, aminus 
     893!CT            CALL prt_ctl_info(charout) 
     894!CT         ENDIF 
     895!CT 
     896!CT         zadefi = MIN( aplus, aminus ) 
     897!CT         IF( zadefi == 0.e0 ) THEN  
     898!CT            erp(:,:) = 0.e0 
     899!CT         ELSE 
     900!CT            erp(:,:) = zadefi * ( zerpplus(:,:) / aplus + zerpminus(:,:) / aminus ) 
     901!CT         ENDIF 
     902!CT 
     903!CT      END IF 
    909904#else 
    910905      ! Rigid-lid (emp=emps=E-P-R+Erp) 
Note: See TracChangeset for help on using the changeset viewer.