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

Ignore:
Timestamp:
2015-02-17T10:06:39+01:00 (9 years ago)
Author:
timgraham
Message:

Merged head of trunk into branch in preparation for putting code back onto the trunk
In working copy ran the command:
svn merge svn+sshtimgraham@…/ipsl/forge/projets/nemo/svn/trunk

Also recompiled NEMO_book.pdf with merged input files

File:
1 edited

Legend:

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

    r4624 r5086  
    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  
    41 #if defined key_lim3 || defined key_cice 
     45   USE sbcwave, ONLY   :  cdn_wave ! wave module 
    4246   USE sbc_ice         ! Surface boundary condition: ice fields 
    43 #endif 
    4447   USE lib_fortran     ! to use key_nosignedzero 
    4548 
     
    5255   PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module 
    5356 
    54    INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read  
     57   INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read 
    5558   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point 
    5659   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point 
     
    6265   INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s) 
    6366   INTEGER , PARAMETER ::   jp_tdif = 9           ! index of tau diff associated to HF tau   (N/m2)   at T-point 
    64     
     67 
    6568   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read) 
    66           
     69 
    6770   !                                             !!! CORE bulk parameters 
    6871   REAL(wp), PARAMETER ::   rhoa =    1.22        ! air density 
     
    7578 
    7679   !                                  !!* Namelist namsbc_core : CORE bulk parameters 
    77    LOGICAL  ::   ln_2m       ! logical flag for height of air temp. and hum 
    7880   LOGICAL  ::   ln_taudif   ! logical flag to use the "mean of stress module - module of mean stress" data 
    7981   REAL(wp) ::   rn_pfac     ! multiplication factor for precipitation 
    8082   REAL(wp) ::   rn_efac     ! multiplication factor for evaporation (clem) 
    8183   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 
    8384   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements 
    8485   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements 
     
    8889#  include "vectopt_loop_substitute.h90" 
    8990   !!---------------------------------------------------------------------- 
    90    !! NEMO/OPA 3.3 , NEMO-consortium (2010)  
     91   !! NEMO/OPA 3.7 , NEMO-consortium (2014) 
    9192   !! $Id$ 
    9293   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    9798      !!--------------------------------------------------------------------- 
    9899      !!                    ***  ROUTINE sbc_blk_core  *** 
    99       !!                    
     100      !! 
    100101      !! ** Purpose :   provide at each time step the surface ocean fluxes 
    101       !!      (momentum, heat, freshwater and runoff)  
     102      !!      (momentum, heat, freshwater and runoff) 
    102103      !! 
    103104      !! ** Method  : (1) READ each fluxes in NetCDF files: 
     
    118119      !! ** Action  :   defined at each time-step at the air-sea interface 
    119120      !!              - utau, vtau  i- and j-component of the wind stress 
    120       !!              - taum, wndm  wind stress and 10m wind modules at T-point 
     121      !!              - taum        wind stress module at T-point 
     122      !!              - wndm        wind speed  module at T-point over free ocean or leads in presence of sea-ice 
    121123      !!              - qns, qsr    non-solar and solar heat fluxes 
    122124      !!              - emp         upward mass flux (evapo. - precip.) 
    123125      !!              - sfx         salt flux due to freezing/melting (non-zero only if ice is present) 
    124126      !!                            (set in limsbc(_2).F90) 
     127      !! 
     128      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008 
     129      !!                   Brodeau et al. Ocean Modelling 2010 
    125130      !!---------------------------------------------------------------------- 
    126131      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    127       !! 
     132      ! 
    128133      INTEGER  ::   ierror   ! return error code 
    129134      INTEGER  ::   ifpr     ! dummy loop indice 
    130135      INTEGER  ::   jfld     ! dummy loop arguments 
    131136      INTEGER  ::   ios      ! Local integer output status for namelist read 
    132       !! 
     137      ! 
    133138      CHARACTER(len=100) ::  cn_dir   !   Root directory for location of core files 
    134139      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read 
     
    136141      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !   "                                 " 
    137142      TYPE(FLD_N) ::   sn_tdif                                 !   "                                 " 
    138       NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac, rn_efac, rn_vfac,  & 
     143      NAMELIST/namsbc_core/ cn_dir , ln_taudif, rn_pfac, rn_efac, rn_vfac,  & 
    139144         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           & 
    140145         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           & 
    141          &                  sn_tdif, rn_zqt , ln_bulk2z, rn_zu 
    142       !!--------------------------------------------------------------------- 
    143  
     146         &                  sn_tdif, rn_zqt, rn_zu 
     147      !!--------------------------------------------------------------------- 
     148      ! 
    144149      !                                         ! ====================== ! 
    145150      IF( kt == nit000 ) THEN                   !  First call kt=nit000  ! 
     
    149154         READ  ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901) 
    150155901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_core in reference namelist', lwp ) 
    151  
     156         ! 
    152157         REWIND( numnam_cfg )              ! Namelist namsbc_core in configuration namelist : CORE bulk parameters 
    153158         READ  ( numnam_cfg, namsbc_core, IOSTAT = ios, ERR = 902 ) 
    154159902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_core in configuration namelist', lwp ) 
    155160 
    156          IF(lwm) WRITE ( numond, namsbc_core ) 
     161         IF(lwm) WRITE( numond, namsbc_core ) 
    157162         !                                         ! 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' )  
     163         IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 )   & 
     164            &   CALL ctl_stop( 'sbc_blk_core: ln_dm2dc can be activated only with daily short-wave forcing' ) 
    160165         IF( ln_dm2dc .AND. sn_qsr%ln_tint ) THEN 
    161166            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' ) 
     167               &         '              ==> We force time interpolation = .false. for qsr' ) 
    163168            sn_qsr%ln_tint = .false. 
    164169         ENDIF 
     
    169174         slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow 
    170175         slf_i(jp_tdif) = sn_tdif 
    171          !                  
     176         ! 
    172177         lhftau = ln_taudif                        ! do we use HF tau information? 
    173178         jfld = jpfld - COUNT( (/.NOT. lhftau/) ) 
     
    191196      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_core( kt, sf, sst_m, ssu_m, ssv_m ) 
    192197 
    193       ! If diurnal cycle is activated, compute a daily mean short waves flux for biogeochemistery  
     198      ! If diurnal cycle is activated, compute a daily mean short waves flux for biogeochemistery 
    194199      IF( ltrcdm2dc )   CALL blk_bio_meanqsr 
    195200 
     
    226231      !!              - qsr     : Solar heat flux over the ocean        (W/m2) 
    227232      !!              - qns     : Non Solar heat flux over the ocean    (W/m2) 
    228       !!              - evap    : Evaporation over the ocean            (kg/m2/s) 
    229233      !!              - emp     : evaporation minus precipitation       (kg/m2/s) 
    230234      !! 
     
    269273      zwnd_j(:,:) = 0.e0 
    270274#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 ! 
     275      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
    275276      DO jj = 2, jpjm1 
    276277         DO ji = fs_2, fs_jpim1   ! vect. opt. 
     
    279280         END DO 
    280281      END DO 
    281 #endif 
    282 #if defined key_vectopt_loop 
    283 !CDIR COLLAPSE 
    284282#endif 
    285283      DO jj = 2, jpjm1 
     
    292290      CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. ) 
    293291      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
    294 !CDIR NOVERRCHK 
    295 !CDIR COLLAPSE 
    296292      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
    297293         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
     
    300296      !      I   Radiative FLUXES                                                     ! 
    301297      ! ----------------------------------------------------------------------------- ! 
    302      
     298 
    303299      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave 
    304300      zztmp = 1. - albo 
     
    306302      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
    307303      ENDIF 
    308 !CDIR COLLAPSE 
    309304      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    310305      ! ----------------------------------------------------------------------------- ! 
     
    313308 
    314309      ! ... specific humidity at SST and IST 
    315 !CDIR NOVERRCHK 
    316 !CDIR COLLAPSE 
    317       zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) )  
     310      zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) ) 
    318311 
    319312      ! ... 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  
     313      CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   & 
     314         &               Cd, Ch, Ce, zt_zu, zq_zu ) 
     315     
    354316      ! ... tau module, i and j component 
    355317      DO jj = 1, jpj 
     
    363325 
    364326      ! ... add the HF tau contribution to the wind stress module? 
    365       IF( lhftau ) THEN  
    366 !CDIR COLLAPSE 
     327      IF( lhftau ) THEN 
    367328         taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 
    368329      ENDIF 
     
    371332      ! ... utau, vtau at U- and V_points, resp. 
    372333      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
     334      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
    373335      DO jj = 1, jpjm1 
    374336         DO ji = 1, fs_jpim1 
    375             utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) 
    376             vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) 
     337            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
     338               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
     339            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) & 
     340               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
    377341         END DO 
    378342      END DO 
     
    380344      CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
    381345 
     346     
    382347      !  Turbulent fluxes over ocean 
    383348      ! ----------------------------- 
    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 
     349      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
     350         !! q_air and t_air are (or "are almost") given at 10m (wind reference height) 
     351         zevap(:,:) = rn_efac*MAX( 0._wp,     rhoa*Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:,1) )*wndm(:,:) ) ! Evaporation 
     352         zqsb (:,:) =                     cpa*rhoa*Ch(:,:)*( zst   (:,:) - sf(jp_tair)%fnow(:,:,1) )*wndm(:,:)   ! Sensible Heat 
    388353      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 
     354         !! q_air and t_air are not given at 10m (wind reference height) 
     355         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
     356         zevap(:,:) = rn_efac*MAX( 0._wp,     rhoa*Ce(:,:)*( zqsatw(:,:) - zq_zu(:,:) )*wndm(:,:) )   ! Evaporation 
     357         zqsb (:,:) =                     cpa*rhoa*Ch(:,:)*( zst   (:,:) - zt_zu(:,:) )*wndm(:,:)     ! Sensible Heat 
    393358      ENDIF 
    394 !CDIR COLLAPSE 
    395359      zqla (:,:) = Lv * zevap(:,:)                                                              ! Latent Heat 
    396360 
     
    409373      !     III    Total FLUXES                                                       ! 
    410374      ! ----------------------------------------------------------------------------- ! 
    411       
    412 !CDIR COLLAPSE 
     375      ! 
    413376      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.) 
    414377         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    415 !CDIR COLLAPSE 
    416378      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar flux 
    417379         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip 
    418380         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
    419381         &     + ( 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                          &    
     382         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          & 
    421383         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
    422          &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic  
     384         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) 
    423385      ! 
    424386      CALL iom_put( "qlw_oce",   zqlw )                 ! output downward longwave heat over the ocean 
     
    442404      ! 
    443405   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   
    498406  
    499407    
     
    518426      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pui      ! ice surface velocity (i- and i- components      [m/s] 
    519427      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid) 
    520       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb     ! ice albedo (clear sky) (alb_ice_cs)               [%] 
     428      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb     ! ice albedo (all skies)                            [%] 
    521429      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_taui   ! i- & j-components of surface ice stress        [N/m2] 
    522430      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid) 
     
    538446      REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2, zcoef_dqlw, zcoef_dqla, zcoef_dqsb 
    539447      REAL(wp) ::   zztmp                                        ! temporary variable 
    540       REAL(wp) ::   zcoef_frca                                   ! fractional cloud amount 
    541448      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f                  ! relative wind module and components at F-point 
    542449      REAL(wp) ::             zwndi_t , zwndj_t                  ! relative wind components at T-point 
     
    562469      zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8) 
    563470      zcoef_dqsb   = rhoa * cpa * Cice 
    564       zcoef_frca   = 1.0  - 0.3 
    565471 
    566472!!gm brutal.... 
     
    579485      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
    580486         !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked) 
    581 !CDIR NOVERRCHK 
    582487         DO jj = 2, jpjm1 
    583488            DO ji = 2, jpim1   ! B grid : NO vector opt 
     
    604509         ! 
    605510      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    606 #if defined key_vectopt_loop 
    607 !CDIR COLLAPSE 
    608 #endif 
    609511         DO jj = 2, jpj 
    610512            DO ji = fs_2, jpi   ! vect. opt. 
     
    614516            END DO 
    615517         END DO 
    616 #if defined key_vectopt_loop 
    617 !CDIR COLLAPSE 
    618 #endif 
    619518         DO jj = 2, jpjm1 
    620519            DO ji = fs_2, fs_jpim1   ! vect. opt. 
     
    635534      DO jl = 1, ijpl                       !  Loop over ice categories  ! 
    636535         !                                  ! ========================== ! 
    637 !CDIR NOVERRCHK 
    638 !CDIR COLLAPSE 
    639536         DO jj = 1 , jpj 
    640 !CDIR NOVERRCHK 
    641537            DO ji = 1, jpi 
    642538               ! ----------------------------! 
     
    648544               p_qsr(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
    649545               ! Long  Wave (lw) 
    650                ! iovino 
    651                IF( ff(ji,jj) .GT. 0._wp ) THEN 
    652                   z_qlw(ji,jj,jl) = ( 0.95 * sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    653                ELSE 
    654                   z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    655                ENDIF 
     546               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    656547               ! lw sensitivity 
    657548               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                                
     
    668559                  &                         * (  11637800. * EXP( -5897.8 / pst(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    669560               ! Latent heat sensitivity for ice (Dqla/Dt) 
    670                p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) ) 
     561               IF( p_qla(ji,jj,jl) > 0._wp ) THEN 
     562                  p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) ) 
     563               ELSE 
     564                  p_dqla(ji,jj,jl) = 0._wp 
     565               ENDIF 
     566 
    671567               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    672568               z_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj) 
     
    676572               ! ----------------------------! 
    677573               ! 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)       
     574               p_qns (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - p_qla (ji,jj,jl) 
    679575               ! 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) )     
     576               p_dqns(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + p_dqla(ji,jj,jl) ) 
    681577            END DO 
    682578            ! 
     
    689585      ! thin surface layer and penetrates inside the ice cover 
    690586      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    691      
    692 !CDIR COLLAPSE 
    693       p_fr1(:,:) = ( 0.18 * ( 1.0 - zcoef_frca ) + 0.35 * zcoef_frca ) 
    694 !CDIR COLLAPSE 
    695       p_fr2(:,:) = ( 0.82 * ( 1.0 - zcoef_frca ) + 0.65 * zcoef_frca ) 
    696         
    697 !CDIR COLLAPSE 
     587      ! 
     588      p_fr1(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     589      p_fr2(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     590      ! 
    698591      p_tpr(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s] 
    699 !CDIR COLLAPSE 
    700592      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  
     593      CALL iom_put( 'snowpre', p_spr * 86400. )                  ! Snow precipitation 
     594      CALL iom_put( 'precip' , p_tpr * 86400. )                  ! Total precipitation 
    703595      ! 
    704596      IF(ln_ctl) THEN 
     
    713605      ENDIF 
    714606 
    715       CALL wrk_dealloc( jpi,jpj, z_wnds_t ) 
    716       CALL wrk_dealloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb )  
     607      CALL wrk_dealloc( jpi,jpj,   z_wnds_t ) 
     608      CALL wrk_dealloc( jpi,jpj,   pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    717609      ! 
    718610      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core') 
    719611      ! 
    720612   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   ) 
     613 
     614 
     615   SUBROUTINE blk_bio_meanqsr 
     616      !!--------------------------------------------------------------------- 
     617      !!                     ***  ROUTINE blk_bio_meanqsr 
     618      !!                      
     619      !! ** Purpose :   provide daily qsr_mean for PISCES when 
     620      !!                analytic diurnal cycle is applied in physic 
     621      !!                 
     622      !! ** Method  :   add part where there is no ice 
     623      !!  
     624      !!--------------------------------------------------------------------- 
     625      IF( nn_timing == 1 )  CALL timing_start('blk_bio_meanqsr') 
     626      ! 
     627      qsr_mean(:,:) = (1. - albo ) *  sf(jp_qsr)%fnow(:,:,1) 
     628      ! 
     629      IF( nn_timing == 1 )  CALL timing_stop('blk_bio_meanqsr') 
     630      ! 
     631   END SUBROUTINE blk_bio_meanqsr 
     632  
     633  
     634   SUBROUTINE blk_ice_meanqsr( palb, p_qsr_mean, pdim ) 
     635      !!--------------------------------------------------------------------- 
     636      !! 
     637      !! ** Purpose :   provide the daily qsr_mean over sea_ice for PISCES when 
     638      !!                analytic diurnal cycle is applied in physic 
     639      !! 
     640      !! ** Method  :   compute qsr 
     641      !!  
     642      !!--------------------------------------------------------------------- 
     643      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb       ! ice albedo (clear sky) (alb_ice_cs)               [%] 
     644      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qsr_mean !     solar heat flux over ice (T-point)         [W/m2] 
     645      INTEGER                   , INTENT(in   ) ::   pdim       ! number of ice categories 
     646      ! 
     647      INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays) 
     648      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     649      REAL(wp) ::   zztmp         ! temporary variable 
     650      !!--------------------------------------------------------------------- 
     651      IF( nn_timing == 1 )  CALL timing_start('blk_ice_meanqsr') 
     652      ! 
     653      ijpl  = pdim                            ! number of ice categories 
     654      zztmp = 1. / ( 1. - albo ) 
     655      !                                     ! ========================== ! 
     656      DO jl = 1, ijpl                       !  Loop over ice categories  ! 
     657         !                                  ! ========================== ! 
     658         DO jj = 1 , jpj 
     659            DO ji = 1, jpi 
     660                  p_qsr_mean(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr_mean(ji,jj) 
     661            END DO 
     662         END DO 
     663      END DO 
     664      ! 
     665      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_meanqsr') 
     666      ! 
     667   END SUBROUTINE blk_ice_meanqsr   
     668 
     669 
     670   SUBROUTINE turb_core_2z( zt, zu, sst, T_zt, q_sat, q_zt, dU,    & 
     671      &                      Cd, Ch, Ce , T_zu, q_zu ) 
    725672      !!---------------------------------------------------------------------- 
    726673      !!                      ***  ROUTINE  turb_core  *** 
    727674      !! 
    728675      !! ** 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 : ??? 
     676      !!                fluxes according to Large & Yeager (2004) and Large & Yeager (2008) 
     677      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
     678      !! 
     679      !! ** Method : Monin Obukhov Similarity Theory  
     680      !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 
     681      !! 
     682      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008 
     683      !! 
     684      !! ** Last update: Laurent Brodeau, June 2014: 
     685      !!    - handles both cases zt=zu and zt/=zu 
     686      !!    - optimized: less 2D arrays allocated and less operations 
     687      !!    - better first guess of stability by checking air-sea difference of virtual temperature 
     688      !!       rather than temperature difference only... 
     689      !!    - added function "cd_neutral_10m" that uses the improved parametrization of  
     690      !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 
     691      !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them 
     692      !!      => 'vkarmn' and 'grav' 
    871693      !!---------------------------------------------------------------------- 
    872694      REAL(wp), INTENT(in   )                     ::   zt       ! height for T_zt and q_zt                   [m] 
     
    876698      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_sat    ! sea surface specific humidity         [kg/kg] 
    877699      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] 
     700      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   dU       ! relative wind module at zu            [m/s] 
    879701      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau) 
    880702      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens) 
     
    882704      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   T_zu     ! air temp. shifted at zu                     [K] 
    883705      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 
     706      ! 
     707      INTEGER ::   j_itt 
     708      INTEGER , PARAMETER ::   nb_itt = 5       ! number of itterations 
     709      LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at different height than U 
     710      ! 
     711      REAL(wp), DIMENSION(:,:), POINTER ::   U_zu          ! relative wind at zu                            [m/s] 
    894712      REAL(wp), DIMENSION(:,:), POINTER ::   Ce_n10        ! 10m neutral latent coefficient 
    895713      REAL(wp), DIMENSION(:,:), POINTER ::   Ch_n10        ! 10m neutral sensible coefficient 
    896714      REAL(wp), DIMENSION(:,:), POINTER ::   sqrt_Cd_n10   ! root square of Cd_n10 
    897715      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] 
    903716      REAL(wp), DIMENSION(:,:), POINTER ::   zeta_u        ! stability parameter at height zu 
    904717      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 
     718      REAL(wp), DIMENSION(:,:), POINTER ::   zpsi_h_u, zpsi_m_u 
     719      REAL(wp), DIMENSION(:,:), POINTER ::   ztmp0, ztmp1, ztmp2 
     720      REAL(wp), DIMENSION(:,:), POINTER ::   stab          ! 1st stability test integer 
    909721      !!---------------------------------------------------------------------- 
    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 
     722 
     723      IF( nn_timing == 1 )  CALL timing_start('turb_core_2z') 
     724     
     725      CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 
     726      CALL wrk_alloc( jpi,jpj, zeta_u, stab ) 
     727      CALL wrk_alloc( jpi,jpj, zpsi_h_u, zpsi_m_u, ztmp0, ztmp1, ztmp2 ) 
     728 
     729      l_zt_equal_zu = .FALSE. 
     730      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     731 
     732      IF( .NOT. l_zt_equal_zu )   CALL wrk_alloc( jpi,jpj, zeta_t ) 
     733 
     734      U_zu = MAX( 0.5 , dU )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
     735 
     736      !! First guess of stability:  
     737      ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt 
     738      stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
     739 
     740      !! Neutral coefficients at 10m: 
     741      IF( ln_cdgw ) THEN      ! wave drag case 
     742         cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) 
     743         ztmp0   (:,:) = cdn_wave(:,:) 
    928744      ELSE 
    929         Cd_n10  = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 )  
     745         ztmp0 = cd_neutral_10m( U_zu ) 
    930746      ENDIF 
    931       sqrt_Cd_n10 = sqrt(Cd_n10) 
     747      sqrt_Cd_n10 = SQRT( ztmp0 ) 
    932748      Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 
    933749      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 
    934  
     750     
    935751      !! 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 
     752      Cd = ztmp0   ;   Ce = Ce_n10   ;   Ch = Ch_n10   ;   sqrt_Cd = sqrt_Cd_n10 
     753 
     754      !! Initializing values at z_u with z_t values: 
     755      T_zu = T_zt   ;   q_zu = q_zt 
    940756 
    941757      !!  * Now starting iteration loop 
    942758      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)) 
     759         ! 
     760         ztmp1 = T_zu - sst   ! Updating air/sea differences 
     761         ztmp2 = q_zu - q_sat  
     762 
     763         ! Updating turbulent scales :   (L&Y 2004 eq. (7)) 
     764         ztmp1  = Ch/sqrt_Cd*ztmp1    ! theta* 
     765         ztmp2  = Ce/sqrt_Cd*ztmp2    ! q* 
     766        
     767         ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu 
     768 
     769         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 
     770         ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu)  
     771         !                                                                     ( Cd*U_zu*U_zu is U*^2 at zu) 
     772 
    951773         !! 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; 
     774         zeta_u   = zu*ztmp0   ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 
     775         zpsi_h_u = psi_h( zeta_u ) 
     776         zpsi_m_u = psi_m( zeta_u ) 
     777        
     778         !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c)) 
     779         IF ( .NOT. l_zt_equal_zu ) THEN 
     780            zeta_t = zt*ztmp0 ;  zeta_t = sign( min(abs(zeta_t),10.0), zeta_t ) 
     781            stab = LOG(zu/zt) - zpsi_h_u + psi_h(zeta_t)  ! stab just used as temp array!!! 
     782            T_zu = T_zt + ztmp1/vkarmn*stab    ! ztmp1 is still theta* 
     783            q_zu = q_zt + ztmp2/vkarmn*stab    ! ztmp2 is still q* 
     784            q_zu = max(0., q_zu) 
     785         END IF 
     786        
     787         IF( ln_cdgw ) THEN      ! surface wave case 
     788            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u )  
     789            Cd      = sqrt_Cd * sqrt_Cd 
    973790         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) 
     791           ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
     792           !   In very rare low-wind conditions, the old way of estimating the 
     793           !   neutral wind speed at 10m leads to a negative value that causes the code 
     794           !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
     795           ztmp0 = MAX( 0.25 , U_zu/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)) ) !  U_n10 
     796           ztmp0 = cd_neutral_10m(ztmp0)                                                 ! Cd_n10 
     797           sqrt_Cd_n10 = sqrt(ztmp0) 
     798        
     799           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                     ! L&Y 2004 eq. (6b) 
     800           stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability 
     801           Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d) 
     802 
     803           !! Update of transfer coefficients: 
     804           ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)   ! L&Y 2004 eq. (10a) 
     805           Cd      = ztmp0 / ( ztmp1*ztmp1 )    
     806           sqrt_Cd = SQRT( Cd ) 
    985807         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          !! 
     808         ! 
     809         ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
     810         ztmp2 = sqrt_Cd / sqrt_Cd_n10 
     811         ztmp1 = 1. + Ch_n10*ztmp0                
     812         Ch  = Ch_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10b) 
     813         ! 
     814         ztmp1 = 1. + Ce_n10*ztmp0                
     815         Ce  = Ce_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
     816         ! 
    996817      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) 
     818 
     819      CALL wrk_dealloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 
     820      CALL wrk_dealloc( jpi,jpj, zeta_u, stab ) 
     821      CALL wrk_dealloc( jpi,jpj, zpsi_h_u, zpsi_m_u, ztmp0, ztmp1, ztmp2 ) 
     822 
     823      IF( .NOT. l_zt_equal_zu ) CALL wrk_dealloc( jpi,jpj, zeta_t ) 
     824 
     825      IF( nn_timing == 1 )  CALL timing_stop('turb_core_2z') 
     826      ! 
     827   END SUBROUTINE turb_core_2z 
     828 
     829 
     830   FUNCTION cd_neutral_10m( zw10 ) 
     831      !!---------------------------------------------------------------------- 
     832      !! Estimate of the neutral drag coefficient at 10m as a function  
     833      !! of neutral wind  speed at 10m 
     834      !! 
     835      !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b) 
     836      !! 
     837      !! Author: L. Brodeau, june 2014 
     838      !!----------------------------------------------------------------------     
     839      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zw10           ! scalar wind speed at 10m (m/s) 
     840      REAL(wp), DIMENSION(jpi,jpj)             ::   cd_neutral_10m 
     841      ! 
     842      REAL(wp), DIMENSION(:,:), POINTER ::   rgt33 
     843      !!----------------------------------------------------------------------     
     844      ! 
     845      CALL wrk_alloc( jpi,jpj, rgt33 ) 
     846      ! 
     847      !! When wind speed > 33 m/s => Cyclone conditions => special treatment 
     848      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1   
     849      cd_neutral_10m = 1.e-3 * ( & 
     850         &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. 
     851         &      + rgt33         *      2.34   )                                                    ! zw10 >= 33. 
     852      ! 
     853      CALL wrk_dealloc( jpi,jpj, rgt33) 
     854      ! 
     855   END FUNCTION cd_neutral_10m 
     856 
     857 
     858   FUNCTION psi_m(pta)   !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    1009859      !------------------------------------------------------------------------------- 
    1010       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 
    1011  
    1012       REAL(wp), PARAMETER :: pi = 3.141592653589793_wp 
     860      ! universal profile stability function for momentum 
     861      !------------------------------------------------------------------------------- 
     862      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta 
     863      ! 
    1013864      REAL(wp), DIMENSION(jpi,jpj)             :: psi_m 
    1014865      REAL(wp), DIMENSION(:,:), POINTER        :: X2, X, stabit 
    1015866      !------------------------------------------------------------------------------- 
    1016  
     867      ! 
    1017868      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  
     869      ! 
     870      X2 = SQRT( ABS( 1. - 16.*pta ) )  ;  X2 = MAX( X2 , 1. )   ;   X = SQRT( X2 ) 
     871      stabit = 0.5 + SIGN( 0.5 , pta ) 
     872      psi_m = -5.*pta*stabit  &                                                          ! Stable 
     873         &    + (1. - stabit)*(2.*LOG((1. + X)*0.5) + LOG((1. + X2)*0.5) - 2.*ATAN(X) + rpi*0.5)  ! Unstable 
     874      ! 
    1024875      CALL wrk_dealloc( jpi,jpj, X2, X, stabit ) 
    1025876      ! 
    1026     END FUNCTION psi_m 
    1027  
    1028  
    1029     FUNCTION psi_h( zta )    !! Psis, L & Y eq. (8c), (8d), (8e) 
     877   END FUNCTION psi_m 
     878 
     879 
     880   FUNCTION psi_h( pta )    !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    1030881      !------------------------------------------------------------------------------- 
    1031       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zta 
     882      ! universal profile stability function for temperature and humidity 
     883      !------------------------------------------------------------------------------- 
     884      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pta 
    1032885      ! 
    1033886      REAL(wp), DIMENSION(jpi,jpj)             ::   psi_h 
    1034       REAL(wp), DIMENSION(:,:), POINTER        :: X2, X, stabit 
     887      REAL(wp), DIMENSION(:,:), POINTER        ::   X2, X, stabit 
    1035888      !------------------------------------------------------------------------------- 
    1036  
     889      ! 
    1037890      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  
     891      ! 
     892      X2 = SQRT( ABS( 1. - 16.*pta ) )   ;   X2 = MAX( X2 , 1. )   ;   X = SQRT( X2 ) 
     893      stabit = 0.5 + SIGN( 0.5 , pta ) 
     894      psi_h = -5.*pta*stabit   &                                       ! Stable 
     895         &    + (1. - stabit)*(2.*LOG( (1. + X2)*0.5 ))                ! Unstable 
     896      ! 
    1044897      CALL wrk_dealloc( jpi,jpj, X2, X, stabit ) 
    1045898      ! 
    1046     END FUNCTION psi_h 
    1047    
     899   END FUNCTION psi_h 
     900 
    1048901   !!====================================================================== 
    1049902END MODULE sbcblk_core 
Note: See TracChangeset for help on using the changeset viewer.