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 for trunk/NEMO/OPA_SRC – NEMO

Changeset 700 for trunk/NEMO/OPA_SRC


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

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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.