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 4675 for branches/2014/dev_CNRS0_blk_core/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 – NEMO

Ignore:
Timestamp:
2014-06-18T17:46:53+02:00 (10 years ago)
Author:
gm
Message:

#1348 CORE bulk simplification and optimization: changes in sbcblk_core and namelist_ref

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_CNRS0_blk_core/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r4624 r4675  
    55   !!===================================================================== 
    66   !! History :  1.0  !  2004-08  (U. Schweckendiek)  Original code 
    7    !!            2.0  !  2005-04  (L. Brodeau, A.M. Treguier) additions:  
     7   !!            2.0  !  2005-04  (L. Brodeau, A.M. Treguier) additions: 
    88   !!                           -  new bulk routine for efficiency 
    99   !!                           -  WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files !!!! 
    10    !!                           -  file names and file characteristics in namelist  
    11    !!                           -  Implement reading of 6-hourly fields    
    12    !!            3.0  !  2006-06  (G. Madec) sbc rewritting    
    13    !!             -   !  2006-12  (L. Brodeau) Original code for TURB_CORE_2Z 
     10   !!                           -  file names and file characteristics in namelist 
     11   !!                           -  Implement reading of 6-hourly fields 
     12   !!            3.0  !  2006-06  (G. Madec) sbc rewritting 
     13   !!             -   !  2006-12  (L. Brodeau) Original code for turb_core_2z 
    1414   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put 
    1515   !!            3.3  !  2010-10  (S. Masson)  add diurnal cycle 
    1616   !!            3.4  !  2011-11  (C. Harris) Fill arrays required by CICE 
     17   !!            3.7  !  2014-06  (L. Brodeau) simplification and optimization of CORE bulk 
    1718   !!---------------------------------------------------------------------- 
    1819 
    1920   !!---------------------------------------------------------------------- 
    20    !!   sbc_blk_core  : bulk formulation as ocean surface boundary condition 
    21    !!                   (forced mode, CORE bulk formulea) 
    22    !!   blk_oce_core  : ocean: computes momentum, heat and freshwater fluxes 
    23    !!   blk_ice_core  : ice  : computes momentum, heat and freshwater fluxes 
    24    !!   turb_core     : computes the CORE turbulent transfer coefficients  
     21   !!   sbc_blk_core    : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea) 
     22   !!   blk_oce_core    : computes momentum, heat and freshwater fluxes over ocean 
     23   !!   blk_ice_core    : computes momentum, heat and freshwater fluxes over ice 
     24   !!   blk_bio_meanqsr : compute daily mean short wave radiation over the ocean 
     25   !!   blk_ice_meanqsr : compute daily mean short wave radiation over the ice 
     26   !!   turb_core_2z    : Computes turbulent transfert coefficients 
     27   !!   cd_neutral_10m  : Estimate of the neutral drag coefficient at 10m 
     28   !!   psi_m           : universal profile stability function for momentum 
     29   !!   psi_h           : universal profile stability function for temperature and humidity 
    2530   !!---------------------------------------------------------------------- 
    2631   USE oce             ! ocean dynamics and tracers 
     
    3843   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3944   USE prtctl          ! Print control 
    40    USE sbcwave,ONLY :  cdn_wave !wave module  
     45   USE sbcwave, ONLY   :  cdn_wave ! wave module 
    4146#if defined key_lim3 || defined key_cice 
    4247   USE sbc_ice         ! Surface boundary condition: ice fields 
     
    5257   PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module 
    5358 
    54    INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read  
     59   INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read 
    5560   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point 
    5661   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point 
     
    6267   INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s) 
    6368   INTEGER , PARAMETER ::   jp_tdif = 9           ! index of tau diff associated to HF tau   (N/m2)   at T-point 
    64     
     69 
    6570   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read) 
    66           
     71 
    6772   !                                             !!! CORE bulk parameters 
    6873   REAL(wp), PARAMETER ::   rhoa =    1.22        ! air density 
     
    7580 
    7681   !                                  !!* Namelist namsbc_core : CORE bulk parameters 
    77    LOGICAL  ::   ln_2m       ! logical flag for height of air temp. and hum 
    7882   LOGICAL  ::   ln_taudif   ! logical flag to use the "mean of stress module - module of mean stress" data 
    7983   REAL(wp) ::   rn_pfac     ! multiplication factor for precipitation 
    8084   REAL(wp) ::   rn_efac     ! multiplication factor for evaporation (clem) 
    8185   REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 
    82    LOGICAL  ::   ln_bulk2z   ! logical flag for case where z(q,t) and z(u) are specified in the namelist 
    8386   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements 
    8487   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements 
     
    8891#  include "vectopt_loop_substitute.h90" 
    8992   !!---------------------------------------------------------------------- 
    90    !! NEMO/OPA 3.3 , NEMO-consortium (2010)  
     93   !! NEMO/OPA 3.7 , NEMO-consortium (2014) 
    9194   !! $Id$ 
    9295   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    97100      !!--------------------------------------------------------------------- 
    98101      !!                    ***  ROUTINE sbc_blk_core  *** 
    99       !!                    
     102      !! 
    100103      !! ** Purpose :   provide at each time step the surface ocean fluxes 
    101       !!      (momentum, heat, freshwater and runoff)  
     104      !!      (momentum, heat, freshwater and runoff) 
    102105      !! 
    103106      !! ** Method  : (1) READ each fluxes in NetCDF files: 
     
    123126      !!              - sfx         salt flux due to freezing/melting (non-zero only if ice is present) 
    124127      !!                            (set in limsbc(_2).F90) 
     128      !! 
     129      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008 
     130      !!                   Brodeau et al. Ocean Modelling 2010 
    125131      !!---------------------------------------------------------------------- 
    126132      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    127       !! 
     133      ! 
    128134      INTEGER  ::   ierror   ! return error code 
    129135      INTEGER  ::   ifpr     ! dummy loop indice 
    130136      INTEGER  ::   jfld     ! dummy loop arguments 
    131137      INTEGER  ::   ios      ! Local integer output status for namelist read 
    132       !! 
     138      ! 
    133139      CHARACTER(len=100) ::  cn_dir   !   Root directory for location of core files 
    134140      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read 
     
    136142      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !   "                                 " 
    137143      TYPE(FLD_N) ::   sn_tdif                                 !   "                                 " 
    138       NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac, rn_efac, rn_vfac,  & 
     144      NAMELIST/namsbc_core/ cn_dir , ln_taudif, rn_pfac, rn_efac, rn_vfac,  & 
    139145         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           & 
    140146         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           & 
    141          &                  sn_tdif, rn_zqt , ln_bulk2z, rn_zu 
    142       !!--------------------------------------------------------------------- 
    143  
     147         &                  sn_tdif, rn_zqt, rn_zu 
     148      !!--------------------------------------------------------------------- 
     149      ! 
    144150      !                                         ! ====================== ! 
    145151      IF( kt == nit000 ) THEN                   !  First call kt=nit000  ! 
     
    149155         READ  ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901) 
    150156901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_core in reference namelist', lwp ) 
    151  
     157         ! 
    152158         REWIND( numnam_cfg )              ! Namelist namsbc_core in configuration namelist : CORE bulk parameters 
    153159         READ  ( numnam_cfg, namsbc_core, IOSTAT = ios, ERR = 902 ) 
    154160902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_core in configuration namelist', lwp ) 
    155161 
    156          IF(lwm) WRITE ( numond, namsbc_core ) 
     162         WRITE ( numond, namsbc_core ) 
    157163         !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing? 
    158          IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 )   &  
    159             &   CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' )  
     164         IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 )   & 
     165            &   CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' ) 
    160166         IF( ln_dm2dc .AND. sn_qsr%ln_tint ) THEN 
    161167            CALL ctl_warn( 'sbc_blk_core: ln_dm2dc is taking care of the temporal interpolation of daily qsr',   & 
    162                  &         '              ==> We force time interpolation = .false. for qsr' ) 
     168               &         '              ==> We force time interpolation = .false. for qsr' ) 
    163169            sn_qsr%ln_tint = .false. 
    164170         ENDIF 
     
    169175         slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow 
    170176         slf_i(jp_tdif) = sn_tdif 
    171          !                  
     177         ! 
    172178         lhftau = ln_taudif                        ! do we use HF tau information? 
    173179         jfld = jpfld - COUNT( (/.NOT. lhftau/) ) 
     
    191197      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_core( kt, sf, sst_m, ssu_m, ssv_m ) 
    192198 
    193       ! If diurnal cycle is activated, compute a daily mean short waves flux for biogeochemistery  
     199      ! If diurnal cycle is activated, compute a daily mean short waves flux for biogeochemistery 
    194200      IF( ltrcdm2dc )   CALL blk_bio_meanqsr 
    195201 
     
    269275      zwnd_j(:,:) = 0.e0 
    270276#if defined key_cyclone 
    271 # if defined key_vectopt_loop 
    272 !CDIR COLLAPSE 
    273 # endif 
    274       CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add Manu ! 
     277      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
    275278      DO jj = 2, jpjm1 
    276279         DO ji = fs_2, fs_jpim1   ! vect. opt. 
     
    279282         END DO 
    280283      END DO 
    281 #endif 
    282 #if defined key_vectopt_loop 
    283 !CDIR COLLAPSE 
    284284#endif 
    285285      DO jj = 2, jpjm1 
     
    292292      CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. ) 
    293293      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
    294 !CDIR NOVERRCHK 
    295 !CDIR COLLAPSE 
    296294      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
    297295         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
     
    300298      !      I   Radiative FLUXES                                                     ! 
    301299      ! ----------------------------------------------------------------------------- ! 
    302      
     300 
    303301      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave 
    304302      zztmp = 1. - albo 
     
    306304      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
    307305      ENDIF 
    308 !CDIR COLLAPSE 
    309306      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    310307      ! ----------------------------------------------------------------------------- ! 
     
    313310 
    314311      ! ... specific humidity at SST and IST 
    315 !CDIR NOVERRCHK 
    316 !CDIR COLLAPSE 
    317       zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) )  
     312      zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) ) 
    318313 
    319314      ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point : 
    320       IF( ln_2m ) THEN 
    321          !! If air temp. and spec. hum. are given at different height (2m) than wind (10m) : 
    322          CALL TURB_CORE_2Z(2.,10., zst   , sf(jp_tair)%fnow,         & 
    323             &                      zqsatw, sf(jp_humi)%fnow, wndm,   & 
    324             &                      Cd    , Ch              , Ce  ,   & 
    325             &                      zt_zu , zq_zu                   ) 
    326       ELSE IF( ln_bulk2z ) THEN 
    327          !! If the height of the air temp./spec. hum. and wind are to be specified by hand : 
    328          IF( rn_zqt == rn_zu ) THEN 
    329             !! If air temp. and spec. hum. are at the same height as wind : 
    330             CALL TURB_CORE_1Z( rn_zu, zst   , sf(jp_tair)%fnow(:,:,1),       & 
    331                &                      zqsatw, sf(jp_humi)%fnow(:,:,1), wndm, & 
    332                &                      Cd    , Ch                     , Ce  ) 
    333          ELSE 
    334             !! If air temp. and spec. hum. are at a different height to wind : 
    335             CALL TURB_CORE_2Z(rn_zqt, rn_zu , zst   , sf(jp_tair)%fnow,         & 
    336                &                              zqsatw, sf(jp_humi)%fnow, wndm,   & 
    337                &                              Cd    , Ch              , Ce  ,   & 
    338                &                              zt_zu , zq_zu                 ) 
    339          ENDIF 
    340       ELSE 
    341          !! If air temp. and spec. hum. are given at same height than wind (10m) : 
    342 !gm bug?  at the compiling phase, add a copy in temporary arrays...  ==> check perf 
    343 !         CALL TURB_CORE_1Z( 10., zst   (:,:), sf(jp_tair)%fnow(:,:),              & 
    344 !            &                    zqsatw(:,:), sf(jp_humi)%fnow(:,:), wndm(:,:),   & 
    345 !            &                    Cd    (:,:),             Ch  (:,:), Ce  (:,:)  ) 
    346 !gm bug 
    347 ! ARPDBG - this won't compile with gfortran. Fix but check performance 
    348 ! as per comment above. 
    349          CALL TURB_CORE_1Z( 10., zst   , sf(jp_tair)%fnow(:,:,1),       & 
    350             &                    zqsatw, sf(jp_humi)%fnow(:,:,1), wndm, & 
    351             &                    Cd    , Ch              , Ce    ) 
    352       ENDIF 
    353  
     315      CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   & 
     316         &               Cd, Ch, Ce, zt_zu, zq_zu ) 
     317     
    354318      ! ... tau module, i and j component 
    355319      DO jj = 1, jpj 
     
    363327 
    364328      ! ... add the HF tau contribution to the wind stress module? 
    365       IF( lhftau ) THEN  
    366 !CDIR COLLAPSE 
     329      IF( lhftau ) THEN 
    367330         taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 
    368331      ENDIF 
     
    380343      CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
    381344 
     345     
    382346      !  Turbulent fluxes over ocean 
    383347      ! ----------------------------- 
    384       IF( ln_2m .OR. ( ln_bulk2z .AND. rn_zqt /= rn_zu ) ) THEN 
    385          ! Values of temp. and hum. adjusted to height of wind must be used 
    386          zevap(:,:) = rn_efac * MAX( 0.e0, rhoa    *Ce(:,:)*( zqsatw(:,:) - zq_zu(:,:) ) * wndm(:,:) )  ! Evaporation 
    387          zqsb (:,:) =                      rhoa*cpa*Ch(:,:)*( zst   (:,:) - zt_zu(:,:) ) * wndm(:,:)     ! Sensible Heat 
     348      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
     349         !! q_air and t_air are (or "are almost") given at 10m (wind reference height) 
     350         zevap(:,:) = rn_efac*MAX( 0._wp,     rhoa*Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:,1) )*wndm(:,:) ) ! Evaporation 
     351         zqsb (:,:) =                     cpa*rhoa*Ch(:,:)*( zst   (:,:) - sf(jp_tair)%fnow(:,:,1) )*wndm(:,:)   ! Sensible Heat 
    388352      ELSE 
    389 !CDIR COLLAPSE 
    390          zevap(:,:) = rn_efac * MAX( 0.e0, rhoa    *Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:,1) ) * wndm(:,:) )   ! Evaporation 
    391 !CDIR COLLAPSE 
    392          zqsb (:,:) =            rhoa*cpa*Ch(:,:)*( zst   (:,:) - sf(jp_tair)%fnow(:,:,1) ) * wndm(:,:)     ! Sensible Heat 
     353         !! q_air and t_air are not given at 10m (wind reference height) 
     354         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
     355         zevap(:,:) = rn_efac*MAX( 0.e0,     rhoa*Ce(:,:)*( zqsatw(:,:) - zq_zu(:,:) )*wndm(:,:) )   ! Evaporation 
     356         zqsb (:,:) =                    cpa*rhoa*Ch(:,:)*( zst   (:,:) - zt_zu(:,:) )*wndm(:,:)     ! Sensible Heat 
    393357      ENDIF 
    394 !CDIR COLLAPSE 
    395358      zqla (:,:) = Lv * zevap(:,:)                                                              ! Latent Heat 
    396359 
     
    409372      !     III    Total FLUXES                                                       ! 
    410373      ! ----------------------------------------------------------------------------- ! 
    411       
    412 !CDIR COLLAPSE 
     374      ! 
    413375      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.) 
    414376         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    415 !CDIR COLLAPSE 
    416377      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar flux 
    417378         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip 
    418379         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
    419380         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair 
    420          &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &    
     381         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          & 
    421382         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
    422383         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic  
     
    442403      ! 
    443404   END SUBROUTINE blk_oce_core 
    444    
    445    SUBROUTINE blk_bio_meanqsr 
    446       !!--------------------------------------------------------------------- 
    447       !!                     ***  ROUTINE blk_bio_meanqsr 
    448       !!                      
    449       !! ** Purpose :   provide daily qsr_mean for PISCES when 
    450       !!                analytic diurnal cycle is applied in physic 
    451       !!                 
    452       !! ** Method  :   add part where there is no ice 
    453       !!  
    454       !!--------------------------------------------------------------------- 
    455       IF( nn_timing == 1 )  CALL timing_start('blk_bio_meanqsr') 
    456  
    457       qsr_mean(:,:) = (1. - albo ) *  sf(jp_qsr)%fnow(:,:,1) 
    458  
    459       IF( nn_timing == 1 )  CALL timing_stop('blk_bio_meanqsr') 
    460  
    461    END SUBROUTINE blk_bio_meanqsr 
    462   
    463   
    464    SUBROUTINE blk_ice_meanqsr(palb,p_qsr_mean,pdim) 
    465       !!--------------------------------------------------------------------- 
    466       !! 
    467       !! ** Purpose :   provide the daily qsr_mean over sea_ice for PISCES when 
    468       !!                analytic diurnal cycle is applied in physic 
    469       !! 
    470       !! ** Method  :   compute qsr 
    471       !!  
    472       !!--------------------------------------------------------------------- 
    473       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb       ! ice albedo (clear sky) (alb_ice_cs)               [%] 
    474       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qsr_mean !     solar heat flux over ice (T-point)         [W/m2] 
    475       INTEGER                   , INTENT(in   ) ::   pdim       ! number of ice categories 
    476       !! 
    477       INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays) 
    478       INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    479       REAL(wp) ::   zztmp         ! temporary variable 
    480       !!--------------------------------------------------------------------- 
    481       IF( nn_timing == 1 )  CALL timing_start('blk_ice_meanqsr') 
    482       ! 
    483       ijpl  = pdim                            ! number of ice categories 
    484       zztmp = 1. / ( 1. - albo ) 
    485       !                                     ! ========================== ! 
    486       DO jl = 1, ijpl                       !  Loop over ice categories  ! 
    487          !                                  ! ========================== ! 
    488          DO jj = 1 , jpj 
    489             DO ji = 1, jpi 
    490                   p_qsr_mean(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr_mean(ji,jj) 
    491             END DO 
    492          END DO 
    493       END DO 
    494       ! 
    495       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_meanqsr') 
    496       ! 
    497    END SUBROUTINE blk_ice_meanqsr   
    498405  
    499406    
     
    579486      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
    580487         !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked) 
    581 !CDIR NOVERRCHK 
    582488         DO jj = 2, jpjm1 
    583489            DO ji = 2, jpim1   ! B grid : NO vector opt 
     
    604510         ! 
    605511      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    606 #if defined key_vectopt_loop 
    607 !CDIR COLLAPSE 
    608 #endif 
    609512         DO jj = 2, jpj 
    610513            DO ji = fs_2, jpi   ! vect. opt. 
     
    614517            END DO 
    615518         END DO 
    616 #if defined key_vectopt_loop 
    617 !CDIR COLLAPSE 
    618 #endif 
    619519         DO jj = 2, jpjm1 
    620520            DO ji = fs_2, fs_jpim1   ! vect. opt. 
     
    635535      DO jl = 1, ijpl                       !  Loop over ice categories  ! 
    636536         !                                  ! ========================== ! 
    637 !CDIR NOVERRCHK 
    638 !CDIR COLLAPSE 
    639537         DO jj = 1 , jpj 
    640 !CDIR NOVERRCHK 
    641538            DO ji = 1, jpi 
    642539               ! ----------------------------! 
     
    676573               ! ----------------------------! 
    677574               ! Downward Non Solar flux 
    678                p_qns (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - p_qla (ji,jj,jl)       
     575               p_qns (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - p_qla (ji,jj,jl) 
    679576               ! Total non solar heat flux sensitivity for ice 
    680                p_dqns(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + p_dqla(ji,jj,jl) )     
     577               p_dqns(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + p_dqla(ji,jj,jl) ) 
    681578            END DO 
    682579            ! 
     
    689586      ! thin surface layer and penetrates inside the ice cover 
    690587      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    691      
    692 !CDIR COLLAPSE 
     588 
    693589      p_fr1(:,:) = ( 0.18 * ( 1.0 - zcoef_frca ) + 0.35 * zcoef_frca ) 
    694 !CDIR COLLAPSE 
    695590      p_fr2(:,:) = ( 0.82 * ( 1.0 - zcoef_frca ) + 0.65 * zcoef_frca ) 
    696         
    697 !CDIR COLLAPSE 
     591 
    698592      p_tpr(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s] 
    699 !CDIR COLLAPSE 
    700593      p_spr(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s] 
    701       CALL iom_put( 'snowpre', p_spr * 86400. )                  ! Snow precipitation  
    702       CALL iom_put( 'precip', p_tpr * 86400. )                   ! Total precipitation  
     594      CALL iom_put( 'snowpre', p_spr * 86400. )                  ! Snow precipitation 
     595      CALL iom_put( 'precip' , p_tpr * 86400. )                  ! Total precipitation 
    703596      ! 
    704597      IF(ln_ctl) THEN 
     
    714607 
    715608      CALL wrk_dealloc( jpi,jpj, z_wnds_t ) 
    716       CALL wrk_dealloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb )  
     609      CALL wrk_dealloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    717610      ! 
    718611      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core') 
    719612      ! 
    720613   END SUBROUTINE blk_ice_core 
    721    
    722  
    723    SUBROUTINE TURB_CORE_1Z(zu, sst, T_a, q_sat, q_a,   & 
    724       &                        dU , Cd , Ch   , Ce   ) 
     614 
     615 
     616   SUBROUTINE blk_bio_meanqsr 
     617      !!--------------------------------------------------------------------- 
     618      !!                     ***  ROUTINE blk_bio_meanqsr 
     619      !!                      
     620      !! ** Purpose :   provide daily qsr_mean for PISCES when 
     621      !!                analytic diurnal cycle is applied in physic 
     622      !!                 
     623      !! ** Method  :   add part where there is no ice 
     624      !!  
     625      !!--------------------------------------------------------------------- 
     626      IF( nn_timing == 1 )  CALL timing_start('blk_bio_meanqsr') 
     627      ! 
     628      qsr_mean(:,:) = (1. - albo ) *  sf(jp_qsr)%fnow(:,:,1) 
     629      ! 
     630      IF( nn_timing == 1 )  CALL timing_stop('blk_bio_meanqsr') 
     631      ! 
     632   END SUBROUTINE blk_bio_meanqsr 
     633  
     634  
     635   SUBROUTINE blk_ice_meanqsr( palb, p_qsr_mean, pdim ) 
     636      !!--------------------------------------------------------------------- 
     637      !! 
     638      !! ** Purpose :   provide the daily qsr_mean over sea_ice for PISCES when 
     639      !!                analytic diurnal cycle is applied in physic 
     640      !! 
     641      !! ** Method  :   compute qsr 
     642      !!  
     643      !!--------------------------------------------------------------------- 
     644      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb       ! ice albedo (clear sky) (alb_ice_cs)               [%] 
     645      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qsr_mean !     solar heat flux over ice (T-point)         [W/m2] 
     646      INTEGER                   , INTENT(in   ) ::   pdim       ! number of ice categories 
     647      ! 
     648      INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays) 
     649      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     650      REAL(wp) ::   zztmp         ! temporary variable 
     651      !!--------------------------------------------------------------------- 
     652      IF( nn_timing == 1 )  CALL timing_start('blk_ice_meanqsr') 
     653      ! 
     654      ijpl  = pdim                            ! number of ice categories 
     655      zztmp = 1. / ( 1. - albo ) 
     656      !                                     ! ========================== ! 
     657      DO jl = 1, ijpl                       !  Loop over ice categories  ! 
     658         !                                  ! ========================== ! 
     659         DO jj = 1 , jpj 
     660            DO ji = 1, jpi 
     661                  p_qsr_mean(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr_mean(ji,jj) 
     662            END DO 
     663         END DO 
     664      END DO 
     665      ! 
     666      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_meanqsr') 
     667      ! 
     668   END SUBROUTINE blk_ice_meanqsr   
     669 
     670 
     671   SUBROUTINE turb_core_2z( zt, zu, sst, T_zt, q_sat, q_zt, dU,    & 
     672      &                      Cd, Ch, Ce , T_zu, q_zu ) 
    725673      !!---------------------------------------------------------------------- 
    726674      !!                      ***  ROUTINE  turb_core  *** 
    727675      !! 
    728676      !! ** Purpose :   Computes turbulent transfert coefficients of surface 
    729       !!                fluxes according to Large & Yeager (2004) 
    730       !! 
    731       !! ** Method  :   I N E R T I A L   D I S S I P A T I O N   M E T H O D 
    732       !!      Momentum, Latent and sensible heat exchange coefficients 
    733       !!      Caution: this procedure should only be used in cases when air 
    734       !!      temperature (T_air), air specific humidity (q_air) and wind (dU) 
    735       !!      are provided at the same height 'zzu'! 
    736       !! 
    737       !! References :   Large & Yeager, 2004 : ??? 
    738       !!---------------------------------------------------------------------- 
    739       REAL(wp)                , INTENT(in   ) ::   zu      ! altitude of wind measurement       [m] 
    740       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   sst     ! sea surface temperature         [Kelvin] 
    741       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   T_a     ! potential air temperature       [Kelvin] 
    742       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   q_sat   ! sea surface specific humidity   [kg/kg] 
    743       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   q_a     ! specific air humidity           [kg/kg] 
    744       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   dU      ! wind module |U(zu)-U(0)|        [m/s] 
    745       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   Cd      ! transfert coefficient for momentum       (tau) 
    746       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   Ch      ! transfert coefficient for temperature (Q_sens) 
    747       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   Ce      ! transfert coefficient for evaporation  (Q_lat) 
    748       !! 
    749       INTEGER :: j_itt 
    750       INTEGER , PARAMETER ::   nb_itt = 3 
    751       REAL(wp), PARAMETER ::   grav   = 9.8   ! gravity                        
    752       REAL(wp), PARAMETER ::   kappa  = 0.4   ! von Karman s constant 
    753  
    754       REAL(wp), DIMENSION(:,:), POINTER  ::   dU10          ! dU                                   [m/s] 
    755       REAL(wp), DIMENSION(:,:), POINTER  ::   dT            ! air/sea temperature differeence      [K] 
    756       REAL(wp), DIMENSION(:,:), POINTER  ::   dq            ! air/sea humidity difference          [K] 
    757       REAL(wp), DIMENSION(:,:), POINTER  ::   Cd_n10        ! 10m neutral drag coefficient 
    758       REAL(wp), DIMENSION(:,:), POINTER  ::   Ce_n10        ! 10m neutral latent coefficient 
    759       REAL(wp), DIMENSION(:,:), POINTER  ::   Ch_n10        ! 10m neutral sensible coefficient 
    760       REAL(wp), DIMENSION(:,:), POINTER  ::   sqrt_Cd_n10   ! root square of Cd_n10 
    761       REAL(wp), DIMENSION(:,:), POINTER  ::   sqrt_Cd       ! root square of Cd 
    762       REAL(wp), DIMENSION(:,:), POINTER  ::   T_vpot        ! virtual potential temperature        [K] 
    763       REAL(wp), DIMENSION(:,:), POINTER  ::   T_star        ! turbulent scale of tem. fluct. 
    764       REAL(wp), DIMENSION(:,:), POINTER  ::   q_star        ! turbulent humidity of temp. fluct. 
    765       REAL(wp), DIMENSION(:,:), POINTER  ::   U_star        ! turb. scale of velocity fluct. 
    766       REAL(wp), DIMENSION(:,:), POINTER  ::   L             ! Monin-Obukov length                  [m] 
    767       REAL(wp), DIMENSION(:,:), POINTER  ::   zeta          ! stability parameter at height zu 
    768       REAL(wp), DIMENSION(:,:), POINTER  ::   U_n10         ! neutral wind velocity at 10m         [m]    
    769       REAL(wp), DIMENSION(:,:), POINTER  ::   xlogt, xct, zpsi_h, zpsi_m 
    770        
    771       INTEGER , DIMENSION(:,:), POINTER  ::   stab          ! 1st guess stability test integer 
    772       !!---------------------------------------------------------------------- 
    773       ! 
    774       IF( nn_timing == 1 )  CALL timing_start('TURB_CORE_1Z') 
    775       ! 
    776       CALL wrk_alloc( jpi,jpj, stab )   ! integer 
    777       CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
    778       CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 
    779  
    780       !! * Start 
    781       !! Air/sea differences 
    782       dU10 = max(0.5, dU)     ! we don't want to fall under 0.5 m/s 
    783       dT = T_a - sst          ! assuming that T_a is allready the potential temp. at zzu 
    784       dq = q_a - q_sat 
    785       !!     
    786       !! Virtual potential temperature 
    787       T_vpot = T_a*(1. + 0.608*q_a) 
    788       !! 
    789       !! Neutral Drag Coefficient 
    790       stab    = 0.5 + sign(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0  
    791       IF  ( ln_cdgw ) THEN 
    792         cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 
    793         Cd_n10(:,:) =   cdn_wave 
    794       ELSE 
    795         Cd_n10  = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a) 
    796       ENDIF 
    797       sqrt_Cd_n10 = sqrt(Cd_n10) 
    798       Ce_n10  = 1.e-3 * ( 34.6 * sqrt_Cd_n10 )               !   L & Y eq. (6b) 
    799       Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) !   L & Y eq. (6c), (6d) 
    800       !! 
    801       !! Initializing transfert coefficients with their first guess neutral equivalents : 
    802       Cd = Cd_n10 ;  Ce = Ce_n10 ;  Ch = Ch_n10 ;  sqrt_Cd = sqrt(Cd) 
    803  
    804       !! * Now starting iteration loop 
    805       DO j_itt=1, nb_itt 
    806          !! Turbulent scales : 
    807          U_star  = sqrt_Cd*dU10                !   L & Y eq. (7a) 
    808          T_star  = Ch/sqrt_Cd*dT               !   L & Y eq. (7b) 
    809          q_star  = Ce/sqrt_Cd*dq               !   L & Y eq. (7c) 
    810  
    811          !! Estimate the Monin-Obukov length : 
    812          L  = (U_star**2)/( kappa*grav*(T_star/T_vpot + q_star/(q_a + 1./0.608)) ) 
    813  
    814          !! Stability parameters : 
    815          zeta  = zu/L ;  zeta   = sign( min(abs(zeta),10.0), zeta ) 
    816          zpsi_h  = psi_h(zeta) 
    817          zpsi_m  = psi_m(zeta) 
    818  
    819          IF  ( ln_cdgw ) THEN 
    820            sqrt_Cd=kappa/((kappa/sqrt_Cd_n10) - zpsi_m) ; Cd=sqrt_Cd*sqrt_Cd; 
    821          ELSE 
    822            !! Shifting the wind speed to 10m and neutral stability : 
    823            U_n10 = dU10*1./(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) !  L & Y eq. (9a) 
    824  
    825            !! Updating the neutral 10m transfer coefficients : 
    826            Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a) 
    827            sqrt_Cd_n10 = sqrt(Cd_n10) 
    828            Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b) 
    829            stab    = 0.5 + sign(0.5,zeta) 
    830            Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab))           !  L & Y eq. (6c), (6d) 
    831  
    832            !! Shifting the neutral  10m transfer coefficients to ( zu , zeta ) : 
    833            !! 
    834            xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10) - zpsi_m) 
    835            Cd  = Cd_n10/(xct*xct) ;  sqrt_Cd = sqrt(Cd) 
    836          ENDIF 
    837          !! 
    838          xlogt = log(zu/10.) - zpsi_h 
    839          !! 
    840          xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10 
    841          Ch  = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct 
    842          !! 
    843          xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10 
    844          Ce  = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct 
    845          !! 
    846       END DO 
    847       !! 
    848       CALL wrk_dealloc( jpi,jpj, stab )   ! integer 
    849       CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
    850       CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 
    851       ! 
    852       IF( nn_timing == 1 )  CALL timing_stop('TURB_CORE_1Z') 
    853       ! 
    854     END SUBROUTINE TURB_CORE_1Z 
    855  
    856  
    857     SUBROUTINE TURB_CORE_2Z(zt, zu, sst, T_zt, q_sat, q_zt, dU, Cd, Ch, Ce, T_zu, q_zu) 
    858       !!---------------------------------------------------------------------- 
    859       !!                      ***  ROUTINE  turb_core  *** 
    860       !! 
    861       !! ** Purpose :   Computes turbulent transfert coefficients of surface  
    862       !!                fluxes according to Large & Yeager (2004). 
    863       !! 
    864       !! ** Method  :   I N E R T I A L   D I S S I P A T I O N   M E T H O D 
    865       !!      Momentum, Latent and sensible heat exchange coefficients 
    866       !!      Caution: this procedure should only be used in cases when air 
    867       !!      temperature (T_air) and air specific humidity (q_air) are at a 
    868       !!      different height to wind (dU). 
    869       !! 
    870       !! References :   Large & Yeager, 2004 : ??? 
     677      !!                fluxes according to Large & Yeager (2004) and Large & Yeager (2008) 
     678      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
     679      !! 
     680      !! ** Method : Monin Obukhov Similarity Theory  
     681      !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 
     682      !! 
     683      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008 
     684      !! 
     685      !! ** Last update: Laurent Brodeau, June 2014: 
     686      !!    - handles both cases zt=zu and zt/=zu 
     687      !!    - optimized: less 2D arrays allocated and less operations 
     688      !!    - better first guess of stability by checking air-sea difference of virtual temperature 
     689      !!       rather than temperature difference only... 
     690      !!    - added function "cd_neutral_10m" that uses the improved parametrization of  
     691      !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 
     692      !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them 
     693      !!      => 'vkarmn' and 'grav' 
    871694      !!---------------------------------------------------------------------- 
    872695      REAL(wp), INTENT(in   )                     ::   zt       ! height for T_zt and q_zt                   [m] 
     
    876699      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_sat    ! sea surface specific humidity         [kg/kg] 
    877700      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                 [kg/kg] 
    878       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   dU       ! relative wind module |U(zu)-U(0)|       [m/s] 
     701      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   dU       ! relative wind module at zu            [m/s] 
    879702      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau) 
    880703      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens) 
     
    882705      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   T_zu     ! air temp. shifted at zu                     [K] 
    883706      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. hum.  shifted at zu               [kg/kg] 
    884  
    885       INTEGER :: j_itt 
    886       INTEGER , PARAMETER :: nb_itt = 5              ! number of itterations 
    887       REAL(wp), PARAMETER ::   grav   = 9.8          ! gravity                        
    888       REAL(wp), PARAMETER ::   kappa  = 0.4          ! von Karman's constant 
    889        
    890       REAL(wp), DIMENSION(:,:), POINTER ::   dU10          ! dU                                [m/s] 
    891       REAL(wp), DIMENSION(:,:), POINTER ::   dT            ! air/sea temperature differeence   [K] 
    892       REAL(wp), DIMENSION(:,:), POINTER ::   dq            ! air/sea humidity difference       [K] 
    893       REAL(wp), DIMENSION(:,:), POINTER ::   Cd_n10        ! 10m neutral drag coefficient 
     707      ! 
     708      INTEGER ::   j_itt 
     709      INTEGER , PARAMETER ::   nb_itt = 5       ! number of itterations 
     710      LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at different height than U 
     711      ! 
     712      REAL(wp), DIMENSION(:,:), POINTER ::   U_zu          ! relative wind at zu                            [m/s] 
    894713      REAL(wp), DIMENSION(:,:), POINTER ::   Ce_n10        ! 10m neutral latent coefficient 
    895714      REAL(wp), DIMENSION(:,:), POINTER ::   Ch_n10        ! 10m neutral sensible coefficient 
    896715      REAL(wp), DIMENSION(:,:), POINTER ::   sqrt_Cd_n10   ! root square of Cd_n10 
    897716      REAL(wp), DIMENSION(:,:), POINTER ::   sqrt_Cd       ! root square of Cd 
    898       REAL(wp), DIMENSION(:,:), POINTER ::   T_vpot        ! virtual potential temperature        [K] 
    899       REAL(wp), DIMENSION(:,:), POINTER ::   T_star        ! turbulent scale of tem. fluct. 
    900       REAL(wp), DIMENSION(:,:), POINTER ::   q_star        ! turbulent humidity of temp. fluct. 
    901       REAL(wp), DIMENSION(:,:), POINTER ::   U_star        ! turb. scale of velocity fluct. 
    902       REAL(wp), DIMENSION(:,:), POINTER ::   L             ! Monin-Obukov length                  [m] 
    903717      REAL(wp), DIMENSION(:,:), POINTER ::   zeta_u        ! stability parameter at height zu 
    904718      REAL(wp), DIMENSION(:,:), POINTER ::   zeta_t        ! stability parameter at height zt 
    905       REAL(wp), DIMENSION(:,:), POINTER ::   U_n10         ! neutral wind velocity at 10m        [m] 
    906       REAL(wp), DIMENSION(:,:), POINTER ::   xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m 
    907  
    908       INTEGER , DIMENSION(:,:), POINTER ::   stab          ! 1st stability test integer 
     719      REAL(wp), DIMENSION(:,:), POINTER ::   zpsi_h_u, zpsi_m_u 
     720      REAL(wp), DIMENSION(:,:), POINTER ::   ztmp0, ztmp1, ztmp2 
     721      REAL(wp), DIMENSION(:,:), POINTER ::   stab          ! 1st stability test integer 
    909722      !!---------------------------------------------------------------------- 
    910       ! 
    911       IF( nn_timing == 1 )  CALL timing_start('TURB_CORE_2Z') 
    912       ! 
    913       CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
    914       CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 
    915       CALL wrk_alloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m ) 
    916       CALL wrk_alloc( jpi,jpj, stab )   ! interger 
    917  
    918       !! Initial air/sea differences 
    919       dU10 = max(0.5, dU)      !  we don't want to fall under 0.5 m/s 
    920       dT = T_zt - sst  
    921       dq = q_zt - q_sat 
    922  
    923       !! Neutral Drag Coefficient : 
    924       stab = 0.5 + sign(0.5,dT)                 ! stab = 1  if dT > 0  -> STABLE 
    925       IF( ln_cdgw ) THEN 
    926         cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 
    927         Cd_n10(:,:) =   cdn_wave 
     723 
     724      IF( nn_timing == 1 )  CALL timing_start('turb_core_2z') 
     725     
     726      CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 
     727      CALL wrk_alloc( jpi,jpj, zeta_u, stab ) 
     728      CALL wrk_alloc( jpi,jpj, zpsi_h_u, zpsi_m_u, ztmp0, ztmp1, ztmp2 ) 
     729 
     730      l_zt_equal_zu = .FALSE. 
     731      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     732 
     733      IF( .NOT. l_zt_equal_zu )   CALL wrk_alloc( jpi,jpj, zeta_t ) 
     734 
     735      U_zu = MAX( 0.5 , dU )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
     736 
     737      !! First guess of stability:  
     738      ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt 
     739      stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
     740 
     741      !! Neutral coefficients at 10m: 
     742      IF( ln_cdgw ) THEN      ! wave drag case 
     743         cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) 
     744         ztmp0   (:,:) = cdn_wave(:,:) 
    928745      ELSE 
    929         Cd_n10  = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 )  
     746         ztmp0 = cd_neutral_10m( U_zu ) 
    930747      ENDIF 
    931       sqrt_Cd_n10 = sqrt(Cd_n10) 
     748      sqrt_Cd_n10 = SQRT( ztmp0 ) 
    932749      Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 
    933750      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 
    934  
     751     
    935752      !! Initializing transf. coeff. with their first guess neutral equivalents : 
    936       Cd = Cd_n10 ;  Ce = Ce_n10 ;  Ch = Ch_n10 ;  sqrt_Cd = sqrt(Cd) 
    937  
    938       !! Initializing z_u values with z_t values : 
    939       T_zu = T_zt ;  q_zu = q_zt 
     753      Cd = ztmp0   ;   Ce = Ce_n10   ;   Ch = Ch_n10   ;   sqrt_Cd = sqrt_Cd_n10 
     754 
     755      !! Initializing values at z_u with z_t values: 
     756      T_zu = T_zt   ;   q_zu = q_zt 
    940757 
    941758      !!  * Now starting iteration loop 
    942759      DO j_itt=1, nb_itt 
    943          dT = T_zu - sst ;  dq = q_zu - q_sat ! Updating air/sea differences 
    944          T_vpot = T_zu*(1. + 0.608*q_zu)    ! Updating virtual potential temperature at zu 
    945          U_star = sqrt_Cd*dU10                ! Updating turbulent scales :   (L & Y eq. (7)) 
    946          T_star  = Ch/sqrt_Cd*dT              ! 
    947          q_star  = Ce/sqrt_Cd*dq              ! 
    948          !! 
    949          L = (U_star*U_star) &                ! Estimate the Monin-Obukov length at height zu 
    950               & / (kappa*grav/T_vpot*(T_star*(1.+0.608*q_zu) + 0.608*T_zu*q_star)) 
     760         ! 
     761         ztmp1 = T_zu - sst   ! Updating air/sea differences 
     762         ztmp2 = q_zu - q_sat  
     763 
     764         ! Updating turbulent scales :   (L&Y 2004 eq. (7)) 
     765         ztmp1  = Ch/sqrt_Cd*ztmp1    ! theta* 
     766         ztmp2  = Ce/sqrt_Cd*ztmp2    ! q* 
     767        
     768         ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu 
     769 
     770         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 
     771         ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu)  
     772         !                                                                     ( Cd*U_zu*U_zu is U*^2 at zu) 
     773 
    951774         !! Stability parameters : 
    952          zeta_u  = zu/L  ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 
    953          zeta_t  = zt/L  ;  zeta_t = sign( min(abs(zeta_t),10.0), zeta_t ) 
    954          zpsi_hu = psi_h(zeta_u) 
    955          zpsi_ht = psi_h(zeta_t) 
    956          zpsi_m  = psi_m(zeta_u) 
    957          !! 
    958          !! Shifting the wind speed to 10m and neutral stability : (L & Y eq.(9a)) 
    959 !        U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u))) 
    960          U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) 
    961          !! 
    962          !! Shifting temperature and humidity at zu :          (L & Y eq. (9b-9c)) 
    963 !        T_zu = T_zt - T_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t)) 
    964          T_zu = T_zt - T_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht) 
    965 !        q_zu = q_zt - q_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t)) 
    966          q_zu = q_zt - q_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht) 
    967          !! 
    968          !! q_zu cannot have a negative value : forcing 0 
    969          stab = 0.5 + sign(0.5,q_zu) ;  q_zu = stab*q_zu 
    970          !! 
    971          IF( ln_cdgw ) THEN 
    972             sqrt_Cd=kappa/((kappa/sqrt_Cd_n10) - zpsi_m) ; Cd=sqrt_Cd*sqrt_Cd; 
     775         zeta_u   = zu*ztmp0   ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 
     776         zpsi_h_u = psi_h( zeta_u ) 
     777         zpsi_m_u = psi_m( zeta_u ) 
     778        
     779         !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c)) 
     780         IF ( .NOT. l_zt_equal_zu ) THEN 
     781            zeta_t = zt*ztmp0 ;  zeta_t = sign( min(abs(zeta_t),10.0), zeta_t ) 
     782            stab = LOG(zu/zt) - zpsi_h_u + psi_h(zeta_t)  ! stab just used as temp array!!! 
     783            T_zu = T_zt + ztmp1/vkarmn*stab    ! ztmp1 is still theta* 
     784            q_zu = q_zt + ztmp2/vkarmn*stab    ! ztmp2 is still q* 
     785            q_zu = max(0., q_zu) 
     786         END IF 
     787        
     788         IF( ln_cdgw ) THEN      ! surface wave case 
     789            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u )  
     790            Cd      = sqrt_Cd * sqrt_Cd 
    973791         ELSE 
    974            !! Updating the neutral 10m transfer coefficients : 
    975            Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)    ! L & Y eq. (6a) 
    976            sqrt_Cd_n10 = sqrt(Cd_n10) 
    977            Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                 ! L & Y eq. (6b) 
    978            stab    = 0.5 + sign(0.5,zeta_u) 
    979            Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) ! L & Y eq. (6c-6d) 
    980            !! 
    981            !! 
    982            !! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) : 
    983            xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)   ! L & Y eq. (10a) 
    984            Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd) 
     792           ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
     793           !   In very rare low-wind conditions, the old way of estimating the 
     794           !   neutral wind speed at 10m leads to a negative value that causes the code 
     795           !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
     796           ztmp0 = MAX( 0.25 , U_zu/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)) ) !  U_n10 
     797           ztmp0 = cd_neutral_10m(ztmp0)                                                 ! Cd_n10 
     798           sqrt_Cd_n10 = sqrt(ztmp0) 
     799        
     800           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                     ! L&Y 2004 eq. (6b) 
     801           stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability 
     802           Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d) 
     803 
     804           !! Update of transfer coefficients: 
     805           ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)   ! L&Y 2004 eq. (10a) 
     806           Cd      = ztmp0 / ( ztmp1*ztmp1 )    
     807           sqrt_Cd = SQRT( Cd ) 
    985808         ENDIF 
    986          !! 
    987          xlogt = log(zu/10.) - zpsi_hu 
    988          !! 
    989          xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10               ! L & Y eq. (10b) 
    990          Ch  = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct 
    991          !! 
    992          xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10               ! L & Y eq. (10c) 
    993          Ce  = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct 
    994          !! 
    995          !! 
     809         ! 
     810         ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
     811         ztmp2 = sqrt_Cd / sqrt_Cd_n10 
     812         ztmp1 = 1. + Ch_n10*ztmp0                
     813         Ch  = Ch_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10b) 
     814         ! 
     815         ztmp1 = 1. + Ce_n10*ztmp0                
     816         Ce  = Ce_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
     817         ! 
    996818      END DO 
    997       !! 
    998       CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
    999       CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 
    1000       CALL wrk_dealloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m ) 
    1001       CALL wrk_dealloc( jpi,jpj, stab )   ! interger 
    1002       ! 
    1003       IF( nn_timing == 1 )  CALL timing_stop('TURB_CORE_2Z') 
    1004       ! 
    1005     END SUBROUTINE TURB_CORE_2Z 
    1006  
    1007  
    1008     FUNCTION psi_m(zta)   !! Psis, L & Y eq. (8c), (8d), (8e) 
     819 
     820      CALL wrk_dealloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 
     821      CALL wrk_dealloc( jpi,jpj, zeta_u, stab ) 
     822      CALL wrk_dealloc( jpi,jpj, zpsi_h_u, zpsi_m_u, ztmp0, ztmp1, ztmp2 ) 
     823 
     824      IF( .NOT. l_zt_equal_zu ) CALL wrk_dealloc( jpi,jpj, zeta_t ) 
     825 
     826      IF( nn_timing == 1 )  CALL timing_stop('turb_core_2z') 
     827      ! 
     828   END SUBROUTINE turb_core_2z 
     829 
     830 
     831   FUNCTION cd_neutral_10m( zw10 ) 
     832      !!---------------------------------------------------------------------- 
     833      !! Estimate of the neutral drag coefficient at 10m as a function  
     834      !! of neutral wind  speed at 10m 
     835      !! 
     836      !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b) 
     837      !! 
     838      !! Author: L. Brodeau, june 2014 
     839      !!----------------------------------------------------------------------     
     840      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zw10           ! scalar wind speed at 10m (m/s) 
     841      REAL(wp), DIMENSION(jpi,jpj)             ::   cd_neutral_10m 
     842      ! 
     843      REAL(wp), DIMENSION(:,:), POINTER ::   rgt33 
     844      !!----------------------------------------------------------------------     
     845      ! 
     846      CALL wrk_alloc( jpi,jpj, rgt33 ) 
     847      ! 
     848      !! When wind speed > 33 m/s => Cyclone conditions => special treatment 
     849      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1   
     850      cd_neutral_10m = 1.e-3 * ( & 
     851         &       (rgt33 + 1._wp)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. 
     852         &      + rgt33         *      2.34   )                                                    ! zw10 >= 33. 
     853      ! 
     854      CALL wrk_dealloc( jpi,jpj, rgt33) 
     855      ! 
     856   END FUNCTION cd_neutral_10m 
     857 
     858 
     859   FUNCTION psi_m(pta)   !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    1009860      !------------------------------------------------------------------------------- 
    1010       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 
    1011  
    1012       REAL(wp), PARAMETER :: pi = 3.141592653589793_wp 
     861      ! universal profile stability function for momentum 
     862      !------------------------------------------------------------------------------- 
     863      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta 
     864      ! 
    1013865      REAL(wp), DIMENSION(jpi,jpj)             :: psi_m 
    1014866      REAL(wp), DIMENSION(:,:), POINTER        :: X2, X, stabit 
    1015867      !------------------------------------------------------------------------------- 
    1016  
     868      ! 
    1017869      CALL wrk_alloc( jpi,jpj, X2, X, stabit ) 
    1018  
    1019       X2 = sqrt(abs(1. - 16.*zta))  ;  X2 = max(X2 , 1.0) ;  X  = sqrt(X2) 
    1020       stabit    = 0.5 + sign(0.5,zta) 
    1021       psi_m = -5.*zta*stabit  &                                                          ! Stable 
    1022          &    + (1. - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2)  ! Unstable  
    1023  
     870      ! 
     871      X2 = SQRT( ABS( 1. - 16.*pta ) )  ;  X2 = MAX( X2 , 1. )   ;   X = SQRT( X2 ) 
     872      stabit = 0.5 + SIGN( 0.5 , pta ) 
     873      psi_m = -5.*pta*stabit  &                                                          ! Stable 
     874         &    + (1. - stabit)*(2.*LOG((1. + X)*0.5) + LOG((1. + X2)*0.5) - 2.*ATAN(X) + rpi*0.5)  ! Unstable 
     875      ! 
    1024876      CALL wrk_dealloc( jpi,jpj, X2, X, stabit ) 
    1025877      ! 
    1026     END FUNCTION psi_m 
    1027  
    1028  
    1029     FUNCTION psi_h( zta )    !! Psis, L & Y eq. (8c), (8d), (8e) 
     878   END FUNCTION psi_m 
     879 
     880 
     881   FUNCTION psi_h( pta )    !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    1030882      !------------------------------------------------------------------------------- 
    1031       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zta 
     883      ! universal profile stability function for temperature and humidity 
     884      !------------------------------------------------------------------------------- 
     885      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pta 
    1032886      ! 
    1033887      REAL(wp), DIMENSION(jpi,jpj)             ::   psi_h 
    1034       REAL(wp), DIMENSION(:,:), POINTER        :: X2, X, stabit 
     888      REAL(wp), DIMENSION(:,:), POINTER        ::   X2, X, stabit 
    1035889      !------------------------------------------------------------------------------- 
    1036  
     890      ! 
    1037891      CALL wrk_alloc( jpi,jpj, X2, X, stabit ) 
    1038  
    1039       X2 = sqrt(abs(1. - 16.*zta))  ;  X2 = max(X2 , 1.) ;  X  = sqrt(X2) 
    1040       stabit    = 0.5 + sign(0.5,zta) 
    1041       psi_h = -5.*zta*stabit  &                                       ! Stable 
    1042          &    + (1. - stabit)*(2.*log( (1. + X2)/2. ))                 ! Unstable 
    1043  
     892      ! 
     893      X2 = SQRT( ABS( 1. - 16.*pta ) )   ;   X2 = MAX( X2 , 1. )   ;   X = SQRT( X2 ) 
     894      stabit = 0.5 + SIGN( 0.5 , pta ) 
     895      psi_h = -5.*pta*stabit   &                                       ! Stable 
     896         &    + (1. - stabit)*(2.*LOG( (1. + X2)*0.5 ))                ! Unstable 
     897      ! 
    1044898      CALL wrk_dealloc( jpi,jpj, X2, X, stabit ) 
    1045899      ! 
    1046     END FUNCTION psi_h 
    1047    
     900   END FUNCTION psi_h 
     901 
    1048902   !!====================================================================== 
    1049903END MODULE sbcblk_core 
Note: See TracChangeset for help on using the changeset viewer.