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 11615 for NEMO – NEMO

Changeset 11615 for NEMO


Ignore:
Timestamp:
2019-09-30T15:18:21+02:00 (5 years ago)
Author:
laurent
Message:

LB: new cool-skin and warm-layer parameterizations for ECMWF and COARE3.6, COARE3.0 uses same CSWL param as for COARE3.6.

Location:
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk
Files:
2 added
1 deleted
8 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/DOM/phycst.F90

    r11111 r11615  
    3636   REAL(wp), PUBLIC ::   omega                       !: earth rotation parameter           [s-1] 
    3737   REAL(wp), PUBLIC ::   ra       = 6371229._wp      !: earth radius                       [m] 
    38    REAL(wp), PUBLIC ::   grav     = 9.80665_wp       !: gravity                            [m/s2]    
     38   REAL(wp), PARAMETER, PUBLIC ::   grav     = 9.80665_wp       !: gravity                            [m/s2]    !LB 
    3939   REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin] 
    4040 
    … …  
    8282   REAL(wp), PARAMETER, PUBLIC :: rho0_a  = 1.2_wp      !: Approx. of density of air                       [kg/m^3] 
    8383   REAL(wp), PARAMETER, PUBLIC :: rho0_w  = 1025._wp    !: Density of sea-water  (ECMWF->1025)             [kg/m^3] 
     84   REAL(wp), PARAMETER, PUBLIC :: roadrw = rho0_a/rho0_w !: Density ratio 
    8485   REAL(wp), PARAMETER, PUBLIC :: rCp0_w  = 4190._wp    !: Specific heat capacity of seawater (ECMWF 4190) [J/K/kg] 
    8586   REAL(wp), PARAMETER, PUBLIC :: rnu0_w  = 1.e-6_wp    !: kinetic viscosity of water                      [m^2/s] 
    … …  
    9091   !                                                    !: the small fraction of downwelling shortwave reflected at the 
    9192   !                                                    !: surface (Lind & Katsaros, 1986) 
    92    REAL(wp), PARAMETER, PUBLIC :: rdct_qsat_salt = 0.98_wp  !: reduction factor on specific humidity at saturation (q_sat(T_s)) due to salt   
     93   REAL(wp), PARAMETER, PUBLIC :: rdct_qsat_salt = 0.98_wp  !: reduction factor on specific humidity at saturation (q_sat(T_s)) due to salt 
     94   REAL(wp), PARAMETER, PUBLIC :: rtt0 = 273.16_wp        !: triple point of temperature    [K] 
     95   REAL(wp), PARAMETER, PUBLIC :: rcst_cs = 16._wp*grav*rho0_w*rCp0_w*rnu0_w*rnu0_w*rnu0_w/(rk0_w*rk0_w) !: for cool-skin parameterizations... 
    9396   !#LB. 
    9497 

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbc_oce.F90

    r11266 r11615  
    153153   !!                     Cool-skin/Warm-layer 
    154154   !!---------------------------------------------------------------------- 
    155    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk       !: sea-surface skin temperature (used if ln_skin==.true.)  [K] !LB 
     155   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk       !: sea-surface skin temperature (used if ln_skin_cs==.true. .OR. ln_skin_wl==.true.)  [K] !LB 
    156156 
    157157    

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90

    r11291 r11615  
    4848   USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003) 
    4949   USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013) 
    50    USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31) 
     50   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 45r1) 
    5151   ! 
    5252   USE iom            ! I/O manager library 
    … …  
    8888   LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    8989   LOGICAL  ::   ln_COARE_3p6   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
    90    LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31) 
     90   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 45r1) 
    9191   ! 
    9292   LOGICAL  ::   ln_taudif      ! logical flag to use the "mean of stress module - module of mean stress" data 
    … …  
    108108 
    109109   !LB: 
    110    LOGICAL  ::   ln_skin        ! use the cool-skin / warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB 
     110   LOGICAL  ::   ln_skin_cs     ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB 
     111   LOGICAL  ::   ln_skin_wl     ! use the warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB 
    111112   LOGICAL  ::   ln_humi_sph    ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB 
    112113   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 
    … …  
    125126   INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    126127   INTEGER, PARAMETER ::   np_COARE_3p6 = 3   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
    127    INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 31) 
     128   INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 45r1) 
    128129 
    129130   !! * Substitutions 
    … …  
    183184         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
    184185         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           & 
    185          &                 ln_skin, ln_humi_sph, ln_humi_dpt, ln_humi_rlh    ! cool-skin / warm-layer !LB 
     186         &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh    ! cool-skin / warm-layer !LB 
    186187      !!--------------------------------------------------------------------- 
    187188      ! 
    … …  
    221222      !LB: 
    222223      !                             !** initialization of the cool-skin / warm-layer parametrization 
    223       IF( ln_skin ) THEN 
     224      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    224225         IF ( ln_NCAR ) CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' ) 
    225226         !                       ! allocate array(s) for cool-skin/warm-layer param. 
    … …  
    246247            CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   & 
    247248               &           '              ==> We force time interpolation = .false. for qsr' ) 
    248             sn_qsr%ln_tint = .FALSE. 
     249            sn_qsr%ln_tint = .false. 
    249250         ENDIF 
    250251      ENDIF 
    … …  
    301302         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0 
    302303         WRITE(numout,*) '      "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 
    303          WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF 
     304         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 45r1)            ln_ECMWF     = ', ln_ECMWF 
    304305         WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif 
    305306         WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt 
    … …  
    317318         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)' 
    318319         CASE( np_COARE_3p6 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 
    319          CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 31)' 
     320         CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 45r1)' 
    320321         END SELECT 
    321322         ! 
    322323         WRITE(numout,*) 
    323          WRITE(numout,*) '      use cool-skin/warm-layer parameterization (SSST)     ln_skin     = ', ln_skin !LB 
     324         WRITE(numout,*) '      use cool-skin  parameterization (SSST)  ln_skin_cs  = ', ln_skin_cs !LB 
     325         WRITE(numout,*) '      use warm-layer parameterization (SSST)  ln_skin_wl  = ', ln_skin_wl !LB 
    324326         ! 
    325327         !LB: 
    … …  
    497499      zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 
    498500 
    499       IF ( ln_skin ) THEN 
     501      IF ( ln_skin_cs .OR. ln_skin_wl ) THEN 
    500502         zqsb(:,:) = zst(:,:) !LB: using array zqsb to backup original zst before skin action 
    501503         zqla(:,:) = zsq(:,:) !LB: using array zqla to backup original zsq before skin action 
    … …  
    522524 
    523525 
    524       IF ( ln_skin ) THEN 
     526      IF ( ln_skin_cs .OR. ln_skin_wl ) THEN 
    525527 
    526528         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    527529 
    528530         CASE( np_COARE_3p0 ) 
    529             CALL turb_coare3p0 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0 
    530                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 
     531            CALL turb_coare3p0 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,                            &  ! COARE v3.0 
     532               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    & 
    531533               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), & 
    532534               &                Tsk_b=tsk(:,:) ) 
    533535 
    534536         CASE( np_COARE_3p6 ) 
    535             CALL turb_coare3p6 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.6 
    536                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 
     537            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITH CSWL options!!!, nsec_day=', nsec_day 
     538            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6 
     539               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    & 
    537540               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), & 
    538                &                Tsk_b=tsk(:,:) ) 
     541               &                isecday_utc=nsec_day, plong=glamt(:,:) ) 
    539542 
    540543         CASE( np_ECMWF     ) 
    541             CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF 
    542                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 
    543                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), & 
    544                &                Tsk_b=tsk(:,:) ) 
     544            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL options!!!, nsec_day=', nsec_day 
     545            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF 
     546               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,  & 
     547               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
    545548 
    546549         CASE DEFAULT 
    547             CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin==.true."' ) 
     550            CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin_*==.true."' ) 
    548551         END SELECT 
    549552 
    … …  
    562565         tsk(:,:) = zst(:,:) 
    563566 
    564       ELSE !IF ( ln_skin ) 
     567      ELSE !IF ( ln_skin_cs .OR. ln_skin_wl ) 
    565568 
    566569 
    567570         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    568571            ! 
    569          CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! NCAR-COREv2 
    570             &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     572         CASE( np_NCAR      ) 
     573            CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! NCAR-COREv2 
     574               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    571575 
    572576         CASE( np_COARE_3p0 ) 
    … …  
    575579               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    576580 
    577          CASE( np_COARE_3p6 )   ;   CALL turb_coare3p6( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.6 
    578             &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     581         CASE( np_COARE_3p6 ) 
     582            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p6" WITHOUT CSWL optional arrays!!!' 
     583            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6 
     584               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    579585 
    580586         CASE( np_ECMWF     ) 
    581587            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!' 
    582             CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF 
     588            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF 
    583589               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    584  
     590             
    585591         CASE DEFAULT 
    586592            CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
    587593         END SELECT 
    588594 
    589       END IF ! IF ( ln_skin ) 
     595      END IF ! IF ( ln_skin_cs .OR. ln_skin_wl ) 
    590596 
    591597 
    … …  
    653659      ! ----------------------------------------------------------------------------- ! 
    654660 
    655       !! LB: now after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin==.TRUE.) 
     661      !! LB: now after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 
    656662      zqlw(:,:) = emiss_w * ( sf(jp_qlw)%fnow(:,:,1) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux 
    657663 
    … …  
    702708      CALL iom_put( 'snowpre', sprecip )                 ! Snow 
    703709      CALL iom_put( 'precip' , tprecip )                 ! Total precipitation 
    704       IF ( ln_skin ) THEN 
     710      IF ( ln_skin_cs .OR. ln_skin_wl ) THEN 
    705711         CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius 
    706712         CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference... 

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90

    r11329 r11615  
    22   !!====================================================================== 
    33   !!                   ***  MODULE  sbcblk_algo_coare3p0  *** 
     4   !! 
     5   !!                After Fairall et al, 2003 
    46   !! Computes: 
    57   !!   * bulk transfer coefficients C_D, C_E and C_H 
    … …  
    79   !!   * the effective bulk wind speed at 10m U_blk 
    810   !!   => all these are used in bulk formulas in sbcblk.F90 
    9    !! 
    10    !!    Using the bulk formulation/param. of COARE v3, Fairall et al., 2003 
    1111   !! 
    1212   !!       Routine turb_coare3p0 maintained and developed in AeroBulk 
    … …  
    3838   USE sbc_oce         ! Surface boundary condition: ocean fields 
    3939   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB 
    40    USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 
     40   USE sbcblk_skin_coare3p6  ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 
    4141 
    4242   IMPLICIT NONE 
    4343   PRIVATE 
    4444 
    45    PUBLIC ::   TURB_COARE3P0   ! called by sbcblk.F90 
     45   PUBLIC :: COARE3P0_INIT, TURB_COARE3P0 
    4646 
    4747   !                                              !! COARE own values for given constants: 
    … …  
    5454CONTAINS 
    5555 
    56    SUBROUTINE turb_coare3p0( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, & 
    57       &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
     56 
     57   SUBROUTINE coare3p0_init(l_use_cs, l_use_wl) 
     58      !!--------------------------------------------------------------------- 
     59      !!                  ***  FUNCTION coare3p0_init  *** 
     60      !! 
     61      !! INPUT : 
     62      !! ------- 
     63      !!    * l_use_cs : use the cool-skin parameterization 
     64      !!    * l_use_wl : use the warm-layer parameterization 
     65      !!--------------------------------------------------------------------- 
     66      LOGICAL , INTENT(in) ::   l_use_cs ! use the cool-skin parameterization 
     67      LOGICAL , INTENT(in) ::   l_use_wl ! use the warm-layer parameterization 
     68      INTEGER :: ierr 
     69      !!--------------------------------------------------------------------- 
     70      IF ( l_use_wl ) THEN 
     71         ierr = 0 
     72         PRINT *, ' *** coare3p0_init: WL => allocating Tau_ac, Qnt_ac, and H_wl :', jpi,jpj 
     73         ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), H_wl(jpi,jpj), STAT=ierr ) 
     74         !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of Tau_ac and Qnt_ac failed!' 
     75         Tau_ac(:,:) = 0._wp 
     76         Qnt_ac(:,:)   = 0._wp 
     77         H_wl(:,:)    = H_wl_max 
     78         PRINT *, ' *** Tau_ac , Qnt_ac, and H_wl allocated!' 
     79      END IF 
     80      !! 
     81      IF ( l_use_cs ) THEN 
     82         ierr = 0 
     83         PRINT *, ' *** coare3p0_init: CS => allocating delta_vl :', jpi,jpj 
     84         ALLOCATE ( delta_vl(jpi,jpj), STAT=ierr ) 
     85         !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of delta_vl and Qnt_ac failed!' 
     86         delta_vl(:,:) = 0.001_wp      ! First guess of zdelta [m] 
     87         PRINT *, ' *** delta_vl allocated!' 
     88      END IF 
     89   END SUBROUTINE coare3p0_init 
     90 
     91 
     92 
     93   SUBROUTINE turb_coare3p0( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl,  & 
     94      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                               & 
    5895      &                      Cdn, Chn, Cen,                      & 
    59       &                      Qsw, rad_lw, slp,                   & 
    60       &                      Tsk_b ) 
     96      &                      Qsw, rad_lw, slp, pdT_cs,                                    & ! optionals for cool-skin (and warm-layer) 
     97      &                      isecday_utc, plong, pdT_wl, Hwl )                             ! optionals for warm-layer only 
    6198      !!---------------------------------------------------------------------- 
    6299      !!                      ***  ROUTINE  turb_coare3p0  *** 
    … …  
    74111      !! INPUT : 
    75112      !! ------- 
     113      !!    *  kt   : current time step (starts at 1) 
    76114      !!    *  zt   : height for temperature and spec. hum. of air            [m] 
    77115      !!    *  zu   : height for wind speed (usually 10m)                     [m] 
    78       !!    *  U_zu : scalar wind speed at zu                                 [m/s] 
    79116      !!    *  t_zt : potential air temperature at zt                         [K] 
    80117      !!    *  q_zt : specific humidity of air at zt                          [kg/kg] 
     118      !!    *  U_zu : scalar wind speed at zu                                 [m/s] 
     119      !!    * l_use_cs : use the cool-skin parameterization 
     120      !!    * l_use_wl : use the warm-layer parameterization 
    81121      !! 
    82122      !! INPUT/OUTPUT: 
    … …  
    87127      !! 
    88128      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg] 
    89       !!              -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True) 
    90       !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False) 
    91       !! 
    92       !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!): 
     129      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True) 
     130      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False) 
     131      !! 
     132      !! OPTIONAL INPUT: 
    93133      !! --------------- 
    94134      !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2] 
    95135      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2] 
    96136      !!    *  slp    : sea-level pressure                                    [Pa] 
    97       !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K] 
     137      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K] 
     138      !!    * isecday_utc: 
     139      !!    *  plong  : longitude array                                       [deg.E] 
     140      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K] 
     141      !!    * Hwl     : depth of warm layer                                   [m] 
    98142      !! 
    99143      !! OUTPUT : 
    … …  
    109153      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    110154      !!---------------------------------------------------------------------------------- 
     155      INTEGER,  INTENT(in   )                     ::   kt       ! current time step 
    111156      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
    112157      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m] 
    … …  
    116161      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg] 
    117162      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s] 
     163      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization 
     164      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization 
    118165      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau) 
    119166      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens) 
    … …  
    127174      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2] 
    128175      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa] 
    129       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa] 
     176      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs 
     177      ! 
     178      INTEGER,  INTENT(in   ), OPTIONAL                     ::   isecday_utc ! current UTC time, counted in second since 00h of the current day 
     179      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   plong    !             [deg.E] 
     180      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K] 
     181      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Hwl      !             [m] 
    130182      ! 
    131183      INTEGER :: j_itt 
    … …  
    141193      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt 
    142194      ! 
    143       ! Cool skin: 
    144       LOGICAL :: l_use_skin = .FALSE. 
    145       REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST 
     195      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & 
     196         &                zsst,   &  ! to back up the initial bulk SST 
     197         &                pdTc,   &  ! SST increment "dT" for cool-skin correction           [K] 
     198         &                pdTw,   &  ! SST increment "dT" for warm layer correction          [K] 
     199         &                zHwl       ! depth of warm-layer [m] 
     200 
     201      ! 
     202      LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE. 
     203      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0.F90' 
     204      CHARACTER(len=128) :: cf_tmp 
    146205      !!---------------------------------------------------------------------------------- 
    147206 
    148       ! Cool skin ? 
    149       IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN 
    150          l_use_skin = .TRUE. 
    151       END IF 
    152       IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin 
     207      IF ( kt == 1 ) CALL COARE3P0_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays 
     208 
    153209 
    154210      l_zt_equal_zu = .FALSE. 
    155211      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
    156  
    157212      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) ) 
    158213 
    159       !! Initialization for cool skin: 
    160       zsst   = T_s    ! save the bulk SST 
    161       IF( l_use_skin ) THEN 
    162          ! First guess for skin temperature: 
    163          IF( PRESENT(Tsk_b) ) THEN 
    164             T_s = Tsk_b 
    165          ELSE 
    166             T_s = T_s - 0.25     ! sst - 0.25 
    167          END IF 
    168          q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 
    169       END IF 
     214      !! Initializations for cool skin and warm layer: 
     215      IF ( l_use_cs ) THEN 
     216         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 
     217            PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!' 
     218            STOP 
     219         END IF 
     220         ALLOCATE ( pdTc(jpi,jpj) ) 
     221         pdTc(:,:) = -0.25_wp  ! First guess of skin correction 
     222      END IF 
     223 
     224      IF ( l_use_wl ) THEN 
     225         IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) .AND. PRESENT(isecday_utc) .AND. PRESENT(plong))) THEN 
     226            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw, slp, isecday_utc & plong to use warm-layer param!' 
     227            STOP 
     228         END IF 
     229         ALLOCATE ( pdTw(jpi,jpj) ) 
     230         IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) ) 
     231      END IF 
     232 
     233      IF ( l_use_cs .OR. l_use_wl ) THEN 
     234         ALLOCATE ( zsst(jpi,jpj) ) 
     235         zsst = T_s ! backing up the bulk SST 
     236         IF( l_use_cs ) T_s = T_s - 0.25   ! First guess of correction 
     237         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!! 
     238      END IF 
     239 
    170240 
    171241      !! First guess of temperature and humidity at height zu: 
    … …  
    190260      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
    191261 
    192       ztmp2  = vkarmn/ztmp0 
    193       Cd     = ztmp2*ztmp2    ! first guess of Cd 
     262      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd 
    194263 
    195264      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 
    … …  
    234303         ztmp1 = u_star*u_star   ! u*^2 
    235304 
    236          !! Update wind at zu taking into acount convection-related wind gustiness: 
    237          ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8): 
    238          ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.) ! square of wind gustiness contribution, ztmp2 == Ug^2 
    239          !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600. 
    240          U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed 
     305         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Fairall et al. 2003, Eq.8): 
     306         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2 
     307         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 
     308         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed 
    241309         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 
    242310 
    … …  
    251319         !! Adjustment the wind at 10m (not needed in the current algo form): 
    252320         !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u)) 
    253           
     321 
    254322         !! Roughness lengthes z0, z0t (z0q = z0t) : 
    255323         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m 
    256324         z0    = alfa_charn_3p0(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6) 
    257          ztmp1 = z0*u_star/znu_a                                           ! Re_r: roughness Reynolds number 
    258          z0t   = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!! 
     325         ztmp1 = ( znu_a / (z0*u_star) )**0.6_wp     ! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific! 
     326         z0t   = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 ) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!! 
    259327 
    260328         !! Turbulent scales at zu : 
    … …  
    267335 
    268336         IF( .NOT. l_zt_equal_zu ) THEN 
    269             ! What's need to be done if zt /= zu 
    270             !! Re-updating temperature and humidity at zu : 
    271             ztmp2 = ztmp0 - psi_h_coare(zeta_t) 
    272             ztmp1 = log(zt/zu) + ztmp2 
     337            !! Re-updating temperature and humidity at zu if zt /= zu : 
     338            ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_coare(zeta_t) 
    273339            t_zu = t_zt - t_star/vkarmn*ztmp1 
    274340            q_zu = q_zt - q_star/vkarmn*ztmp1 
    275341         END IF 
    276342 
    277          !! SKIN related part 
    278          !! ----------------- 
    279          IF( l_use_skin ) THEN 
    280             !! compute transfer coefficients at zu : lolo: verifier... 
    281             ztmp0 = u_star/U_blk 
    282             Ch   = ztmp0*t_star/dt_zu 
    283             Ce   = ztmp0*q_star/dq_zu 
    284             ! Non-Solar heat flux to the ocean: 
    285             ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10 
    286             ztmp2 = T_s*T_s 
    287             ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux 
    288                &       +      emiss_w*(rad_lw - stefan*ztmp2*ztmp2)                         ! Net longwave flux 
    289             !!         => "ztmp1" is the net non-solar surface heat flux ! 
    290             !! Updating the values of the skin temperature T_s and q_s : 
    291             CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!) 
    292             q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp)  ! 200 -> just to avoid numerics problem on masked regions if silly values are given 
    293          END IF 
    294  
    295          IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN 
     343 
     344         IF( l_use_cs ) THEN 
     345            !! Cool-skin contribution 
     346 
     347            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 
     348               &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
     349 
     350            CALL CS_COARE3P6( Qsw, ztmp1, u_star, zsst, ztmp2,  pdTc )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2 
     351 
     352            T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1) 
     353            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)*tmask(:,:,1) 
     354            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
     355 
     356         END IF 
     357 
     358         IF( l_use_wl ) THEN 
     359            !! Warm-layer contribution 
     360            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 
     361               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
     362            !! In WL_COARE3P6 or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 
     363            IF (PRESENT(Hwl)) THEN 
     364               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl ) 
     365            ELSE 
     366               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw ) 
     367            END IF 
     368            !! Updating T_s and q_s !!! 
     369            T_s(:,:) = zsst(:,:) + pdTw(:,:)*tmask(:,:,1) 
     370            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)*tmask(:,:,1) 
     371            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
     372            !LOLO: 
     373            !IF ( j_itt == nb_itt) THEN 
     374            !   WRITE(cf_tmp,'("Qnt_ac_k",i5.5,"_p",i4.4,".nc")') kt, narea 
     375            !   CALL DUMP_FIELD(REAL( Qnt_ac*tmask(:,:,1) , 4), TRIM(cf_tmp), 'Qnt_ac') 
     376            !   WRITE(cf_tmp,  '("pdTw_k",i5.5,"_p",i4.4,".nc")') kt, narea 
     377            !   CALL DUMP_FIELD(REAL( pdTw*tmask(:,:,1) , 4), TRIM(cf_tmp), 'pdTw') 
     378            !END IF 
     379            !LOLO. 
     380         END IF 
     381 
     382 
     383         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 
    296384            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
    297385            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    … …  
    312400 
    313401      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 
     402 
     403      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc*tmask(:,:,1) 
     404      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw*tmask(:,:,1) 
     405 
     406      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
     407      IF (          l_use_cs      ) DEALLOCATE ( pdTc ) 
     408      IF (          l_use_wl      ) THEN 
     409         DEALLOCATE ( pdTw ) 
     410         IF (PRESENT(Hwl)) DEALLOCATE ( zHwl ) 
     411      END IF 
    314412 
    315413   END SUBROUTINE turb_coare3p0 

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90

    r11329 r11615  
    3838   USE sbc_oce         ! Surface boundary condition: ocean fields 
    3939   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB 
    40    USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 
     40   USE sbcblk_skin_coare3p6  ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 
    4141 
    4242   IMPLICIT NONE 
    4343   PRIVATE 
    4444 
    45    PUBLIC ::   TURB_COARE3P6   ! called by sbcblk.F90 
     45   PUBLIC :: COARE3P6_INIT, TURB_COARE3P6 
    4646 
    4747   !                                              !! COARE own values for given constants: 
    … …  
    5454CONTAINS 
    5555 
    56    SUBROUTINE turb_coare3p6( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, & 
    57       &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
     56 
     57   SUBROUTINE coare3p6_init(l_use_cs, l_use_wl) 
     58      !!--------------------------------------------------------------------- 
     59      !!                  ***  FUNCTION coare3p6_init  *** 
     60      !! 
     61      !! INPUT : 
     62      !! ------- 
     63      !!    * l_use_cs : use the cool-skin parameterization 
     64      !!    * l_use_wl : use the warm-layer parameterization 
     65      !!--------------------------------------------------------------------- 
     66      LOGICAL , INTENT(in) ::   l_use_cs ! use the cool-skin parameterization 
     67      LOGICAL , INTENT(in) ::   l_use_wl ! use the warm-layer parameterization 
     68      INTEGER :: ierr 
     69      !!--------------------------------------------------------------------- 
     70      IF ( l_use_wl ) THEN 
     71         ierr = 0 
     72         PRINT *, ' *** coare3p6_init: WL => allocating Tau_ac, Qnt_ac, and H_wl :', jpi,jpj 
     73         ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), H_wl(jpi,jpj), STAT=ierr ) 
     74         !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of Tau_ac and Qnt_ac failed!' 
     75         Tau_ac(:,:) = 0._wp 
     76         Qnt_ac(:,:)   = 0._wp 
     77         H_wl(:,:)    = H_wl_max 
     78         PRINT *, ' *** Tau_ac , Qnt_ac, and H_wl allocated!' 
     79      END IF 
     80      !! 
     81      IF ( l_use_cs ) THEN 
     82         ierr = 0 
     83         PRINT *, ' *** coare3p6_init: CS => allocating delta_vl :', jpi,jpj 
     84         ALLOCATE ( delta_vl(jpi,jpj), STAT=ierr ) 
     85         !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of delta_vl and Qnt_ac failed!' 
     86         delta_vl(:,:) = 0.001_wp      ! First guess of zdelta [m] 
     87         PRINT *, ' *** delta_vl allocated!' 
     88      END IF 
     89   END SUBROUTINE coare3p6_init 
     90 
     91 
     92 
     93   SUBROUTINE turb_coare3p6( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl,  & 
     94      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                               & 
    5895      &                      Cdn, Chn, Cen,                      & 
    59       &                      Qsw, rad_lw, slp,                   & 
    60       &                      Tsk_b ) 
     96      &                      Qsw, rad_lw, slp, pdT_cs,                                    & ! optionals for cool-skin (and warm-layer) 
     97      &                      isecday_utc, plong, pdT_wl, Hwl )                             ! optionals for warm-layer only 
    6198      !!---------------------------------------------------------------------- 
    6299      !!                      ***  ROUTINE  turb_coare3p6  *** 
    … …  
    74111      !! INPUT : 
    75112      !! ------- 
     113      !!    *  kt   : current time step (starts at 1) 
    76114      !!    *  zt   : height for temperature and spec. hum. of air            [m] 
    77115      !!    *  zu   : height for wind speed (usually 10m)                     [m] 
    78       !!    *  U_zu : scalar wind speed at zu                                 [m/s] 
    79116      !!    *  t_zt : potential air temperature at zt                         [K] 
    80117      !!    *  q_zt : specific humidity of air at zt                          [kg/kg] 
     118      !!    *  U_zu : scalar wind speed at zu                                 [m/s] 
     119      !!    * l_use_cs : use the cool-skin parameterization 
     120      !!    * l_use_wl : use the warm-layer parameterization 
    81121      !! 
    82122      !! INPUT/OUTPUT: 
    … …  
    90130      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False) 
    91131      !! 
    92       !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!): 
     132      !! OPTIONAL INPUT: 
    93133      !! --------------- 
    94134      !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2] 
    95135      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2] 
    96136      !!    *  slp    : sea-level pressure                                    [Pa] 
    97       !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K] 
     137      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K] 
     138      !!    * isecday_utc: 
     139      !!    *  plong  : longitude array                                       [deg.E] 
     140      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K] 
     141      !!    * Hwl     : depth of warm layer                                   [m] 
    98142      !! 
    99143      !! OUTPUT : 
    … …  
    109153      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    110154      !!---------------------------------------------------------------------------------- 
     155      INTEGER,  INTENT(in   )                     ::   kt       ! current time step 
    111156      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
    112157      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m] 
    … …  
    116161      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg] 
    117162      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s] 
     163      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization 
     164      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization 
    118165      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau) 
    119166      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens) 
    … …  
    127174      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2] 
    128175      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa] 
    129       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa] 
     176      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs 
     177      ! 
     178      INTEGER,  INTENT(in   ), OPTIONAL                     ::   isecday_utc ! current UTC time, counted in second since 00h of the current day 
     179      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   plong    !             [deg.E] 
     180      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K] 
     181      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Hwl      !             [m] 
    130182      ! 
    131183      INTEGER :: j_itt 
    … …  
    141193      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt 
    142194      ! 
    143       ! Cool skin: 
    144       LOGICAL :: l_use_skin = .FALSE. 
    145       REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST 
     195      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & 
     196         &                zsst,   &  ! to back up the initial bulk SST 
     197         &                pdTc,   &  ! SST increment "dT" for cool-skin correction           [K] 
     198         &                pdTw,   &  ! SST increment "dT" for warm layer correction          [K] 
     199         &                zHwl       ! depth of warm-layer [m] 
     200 
     201      ! 
     202      LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE. 
     203      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p6.F90' 
     204      CHARACTER(len=128) :: cf_tmp 
    146205      !!---------------------------------------------------------------------------------- 
    147206 
    148       ! Cool skin ? 
    149       IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN 
    150          l_use_skin = .TRUE. 
    151       END IF 
    152       IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin 
     207      IF ( kt == 1 ) CALL COARE3P6_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays 
     208 
    153209 
    154210      l_zt_equal_zu = .FALSE. 
    155211      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
    156  
    157212      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) ) 
    158213 
    159       !! Initialization for cool skin: 
    160       zsst   = T_s    ! save the bulk SST 
    161       IF( l_use_skin ) THEN 
    162          ! First guess for skin temperature: 
    163          IF( PRESENT(Tsk_b) ) THEN 
    164             T_s = Tsk_b 
    165          ELSE 
    166             T_s = T_s - 0.25     ! sst - 0.25 
    167          END IF 
    168          q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 
    169       END IF 
     214      !! Initializations for cool skin and warm layer: 
     215      IF ( l_use_cs ) THEN 
     216         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 
     217            PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!' 
     218            STOP 
     219         END IF 
     220         ALLOCATE ( pdTc(jpi,jpj) ) 
     221         pdTc(:,:) = -0.25_wp  ! First guess of skin correction 
     222      END IF 
     223 
     224      IF ( l_use_wl ) THEN 
     225         IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) .AND. PRESENT(isecday_utc) .AND. PRESENT(plong))) THEN 
     226            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw, slp, isecday_utc & plong to use warm-layer param!' 
     227            STOP 
     228         END IF 
     229         ALLOCATE ( pdTw(jpi,jpj) ) 
     230         IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) ) 
     231      END IF 
     232 
     233      IF ( l_use_cs .OR. l_use_wl ) THEN 
     234         ALLOCATE ( zsst(jpi,jpj) ) 
     235         zsst = T_s ! backing up the bulk SST 
     236         IF( l_use_cs ) T_s = T_s - 0.25   ! First guess of correction 
     237         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!! 
     238      END IF 
     239 
    170240 
    171241      !! First guess of temperature and humidity at height zu: 
    … …  
    190260      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
    191261 
    192       ztmp2  = vkarmn/ztmp0 
    193       Cd     = ztmp2*ztmp2    ! first guess of Cd 
     262      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd 
    194263 
    195264      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 
    … …  
    234303         ztmp1 = u_star*u_star   ! u*^2 
    235304 
    236          !! Update wind at zu taking into acount convection-related wind gustiness: 
    237          ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8): 
    238          ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.) ! square of wind gustiness contribution, ztmp2 == Ug^2 
    239          !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600. 
    240          U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed 
     305         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Fairall et al. 2003, Eq.8): 
     306         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2 
     307         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 
     308         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed 
    241309         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 
    242310 
    … …  
    251319         !! Adjustment the wind at 10m (not needed in the current algo form): 
    252320         !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u)) 
    253           
     321 
    254322         !! Roughness lengthes z0, z0t (z0q = z0t) : 
    255323         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m 
    256          z0    = alfa_charn_3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6) 
    257          ztmp1 = z0*u_star/znu_a                                           ! Re_r: roughness Reynolds number 
    258          z0t   = MIN( 1.6E-4_wp , 5.8E-5_wp*ztmp1**(-0.72_wp))   ! COARE 3.6 
     324         z0    = alfa_charn_3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6) [ ztmp1==u*^2 ] 
     325         ztmp1 = ( znu_a / (z0*u_star) )**0.72_wp     ! COARE3.6-specific! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific! 
     326         z0t   = MIN( 1.6E-4_wp , 5.8E-5_wp*ztmp1 )   ! COARE3.6-specific! 
    259327 
    260328         !! Turbulent scales at zu : 
    … …  
    267335 
    268336         IF( .NOT. l_zt_equal_zu ) THEN 
    269             ! What's need to be done if zt /= zu 
    270             !! Re-updating temperature and humidity at zu : 
    271             ztmp2 = ztmp0 - psi_h_coare(zeta_t) 
    272             ztmp1 = log(zt/zu) + ztmp2 
     337            !! Re-updating temperature and humidity at zu if zt /= zu : 
     338            ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_coare(zeta_t) 
    273339            t_zu = t_zt - t_star/vkarmn*ztmp1 
    274340            q_zu = q_zt - q_star/vkarmn*ztmp1 
    275341         END IF 
    276342 
    277          !! SKIN related part 
    278          !! ----------------- 
    279          IF( l_use_skin ) THEN 
    280             !! compute transfer coefficients at zu : lolo: verifier... 
    281             ztmp0 = u_star/U_blk 
    282             Ch   = ztmp0*t_star/dt_zu 
    283             Ce   = ztmp0*q_star/dq_zu 
    284             ! Non-Solar heat flux to the ocean: 
    285             ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10 
    286             ztmp2 = T_s*T_s 
    287             ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux 
    288                &       +      emiss_w*(rad_lw - stefan*ztmp2*ztmp2)                         ! Net longwave flux 
    289             !!         => "ztmp1" is the net non-solar surface heat flux ! 
    290             !! Updating the values of the skin temperature T_s and q_s : 
    291             CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!) 
    292             q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp)  ! 200 -> just to avoid numerics problem on masked regions if silly values are given 
    293          END IF 
    294  
    295          IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN 
     343 
     344         IF( l_use_cs ) THEN 
     345            !! Cool-skin contribution 
     346 
     347            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 
     348               &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
     349 
     350            CALL CS_COARE3P6( Qsw, ztmp1, u_star, zsst, ztmp2,  pdTc )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2 
     351 
     352            T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1) 
     353            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)*tmask(:,:,1) 
     354            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
     355 
     356         END IF 
     357 
     358         IF( l_use_wl ) THEN 
     359            !! Warm-layer contribution 
     360            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 
     361               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
     362            !! In WL_COARE3P6 or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 
     363            IF (PRESENT(Hwl)) THEN 
     364               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl ) 
     365            ELSE 
     366               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw ) 
     367            END IF 
     368            !! Updating T_s and q_s !!! 
     369            T_s(:,:) = zsst(:,:) + pdTw(:,:)*tmask(:,:,1) 
     370            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)*tmask(:,:,1) 
     371            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
     372            !LOLO: 
     373            !IF ( j_itt == nb_itt) THEN 
     374            !   WRITE(cf_tmp,'("Qnt_ac_k",i5.5,"_p",i4.4,".nc")') kt, narea 
     375            !   CALL DUMP_FIELD(REAL( Qnt_ac*tmask(:,:,1) , 4), TRIM(cf_tmp), 'Qnt_ac') 
     376            !   WRITE(cf_tmp,  '("pdTw_k",i5.5,"_p",i4.4,".nc")') kt, narea 
     377            !   CALL DUMP_FIELD(REAL( pdTw*tmask(:,:,1) , 4), TRIM(cf_tmp), 'pdTw') 
     378            !END IF 
     379            !LOLO. 
     380         END IF 
     381 
     382 
     383         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 
    296384            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
    297385            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    … …  
    312400 
    313401      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 
     402 
     403      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc*tmask(:,:,1) 
     404      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw*tmask(:,:,1) 
     405 
     406      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
     407      IF (          l_use_cs      ) DEALLOCATE ( pdTc ) 
     408      IF (          l_use_wl      ) THEN 
     409         DEALLOCATE ( pdTw ) 
     410         IF (PRESENT(Hwl)) DEALLOCATE ( zHwl ) 
     411      END IF 
    314412 
    315413   END SUBROUTINE turb_coare3p6 

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r11329 r11615  
    3939   USE sbc_oce         ! Surface boundary condition: ocean fields 
    4040   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB 
    41    USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 
     41   USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB 
    4242 
    4343   IMPLICIT NONE 
    … …  
    6060CONTAINS 
    6161 
    62    SUBROUTINE turb_ecmwf( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, & 
     62   SUBROUTINE turb_ecmwf( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 
    6363      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
    6464      &                      Cdn, Chn, Cen,                      & 
    65       &                      Qsw, rad_lw, slp,                   & 
    66       &                      Tsk_b ) 
     65      &                      Qsw, rad_lw, slp, pdT_cs,           & ! optionals for cool-skin (and warm-layer) 
     66      &                      pdT_wl )                              ! optionals for warm-layer only 
    6767      !!---------------------------------------------------------------------- 
    68        !!                      ***  ROUTINE  turb_ecmwf  *** 
     68      !!                      ***  ROUTINE  turb_ecmwf  *** 
    6969      !! 
    7070      !! ** Purpose :   Computes turbulent transfert coefficients of surface 
    71       !!                fluxes according to IFS doc. (cycle 40) 
     71      !!                fluxes according to IFS doc. (cycle 45r1) 
    7272      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
    7373      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas 
    … …  
    8282      !!    *  zt   : height for temperature and spec. hum. of air            [m] 
    8383      !!    *  zu   : height for wind speed (usually 10m)                     [m] 
    84       !!    *  U_zu : scalar wind speed at zu                                 [m/s] 
    8584      !!    *  t_zt : potential air temperature at zt                         [K] 
    8685      !!    *  q_zt : specific humidity of air at zt                          [kg/kg] 
     86      !!    *  U_zu : scalar wind speed at zu                                 [m/s] 
     87      !!    * l_use_cs : use the cool-skin parameterization 
     88      !!    * l_use_wl : use the warm-layer parameterization 
    8789      !! 
    8890      !! INPUT/OUTPUT: 
    … …  
    101103      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2] 
    102104      !!    *  slp    : sea-level pressure                                    [Pa] 
    103       !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K] 
     105      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K] 
     106      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K] 
    104107      !! 
    105108      !! OUTPUT : 
    … …  
    122125      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg] 
    123126      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s] 
     127      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization 
     128      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization 
    124129      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau) 
    125130      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens) 
    … …  
    133138      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2] 
    134139      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa] 
    135       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa] 
     140      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs 
     141      ! 
     142      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K] 
    136143      ! 
    137144      INTEGER :: j_itt 
    … …  
    145152         &  z0, z0t, z0q 
    146153      ! 
    147       ! Cool skin: 
    148       LOGICAL :: l_use_skin = .FALSE. 
    149       REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST 
     154      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & 
     155         &                zsst,   &  ! to back up the initial bulk SST 
     156         &                pdTc,   &  ! SST increment "dT" for cool-skin correction           [K] 
     157         &                pdTw       ! SST increment "dT" for warm layer correction          [K] 
    150158      ! 
    151159      REAL(wp), DIMENSION(jpi,jpj) ::   func_m, func_h 
    152160      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2 
    153       !!---------------------------------------------------------------------------------- 
    154  
    155       ! Cool skin ? 
    156       IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN 
    157          l_use_skin = .TRUE. 
    158       END IF 
    159       IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin 
    160  
    161       ! Identical first gess as in COARE, with IFS parameter values though... 
    162       ! 
     161      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90' 
     162      !!---------------------------------------------------------------------------------- 
     163 
    163164      l_zt_equal_zu = .FALSE. 
    164165      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
    165166 
    166       !! Initialization for cool skin: 
    167       zsst   = T_s    ! save the bulk SST 
    168       IF( l_use_skin ) THEN 
    169          ! First guess for skin temperature: 
    170          IF( PRESENT(Tsk_b) ) THEN 
    171             T_s = Tsk_b 
    172          ELSE 
    173             T_s = T_s - 0.25     ! sst - 0.25 
    174          END IF 
    175          q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 
     167      !! Initializations for cool skin and warm layer: 
     168      IF ( l_use_cs ) THEN 
     169         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 
     170            PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!' 
     171            STOP 
     172         END IF 
     173         ALLOCATE ( pdTc(jpi,jpj) ) 
     174         pdTc(:,:) = -0.25_wp  ! First guess of skin correction 
    176175      END IF 
    177176 
     177      IF ( l_use_wl ) THEN 
     178         IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) THEN 
     179            PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!' 
     180            STOP 
     181         END IF 
     182         ALLOCATE ( pdTw(jpi,jpj) ) 
     183         pdTw(:,:) = 0._wp 
     184      END IF 
     185 
     186      IF ( l_use_cs .OR. l_use_wl ) THEN 
     187         ALLOCATE ( zsst(jpi,jpj) ) 
     188         zsst = T_s ! backing up the bulk SST 
     189         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction 
     190         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!! 
     191      END IF 
     192 
     193 
     194      ! Identical first gess as in COARE, with IFS parameter values though... 
     195      ! 
    178196      !! First guess of temperature and humidity at height zu: 
    179197      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions... 
    … …  
    197215      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
    198216 
    199       ztmp2  = vkarmn/ztmp0 
    200       Cd     = ztmp2*ztmp2    ! first guess of Cd 
     217      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd 
    201218 
    202219      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 
    … …  
    265282         z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp) 
    266283 
    267          !! Update wind at 10m taking into acount convection-related wind gustiness: 
    268          !! => Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8 
    269          ! Only true when unstable (L<0) => when ztmp0 < 0 => - !!! 
    270          ztmp2 = ztmp2 * ( MAX(-zi0*Linv/vkarmn , 0._wp))**(2._wp/3._wp) ! => w*^2  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 
    271          !! => equivalent using Beta=1 (gustiness parameter, 1.25 for COARE, also zi0=600 in COARE..) 
    272          U_blk = MAX( SQRT(U_zu*U_zu + ztmp2) , 0.2_wp )              ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1 
     284         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8) 
     285         ztmp2 = Beta0*Beta0*ztmp2*(MAX(-zi0*Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 
     286         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 
     287         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed 
    273288         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 
    274289 
    … …  
    303318         func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 
    304319 
    305          !! SKIN related part 
    306          !! ----------------- 
    307          IF( l_use_skin ) THEN 
    308             !! compute transfer coefficients at zu : lolo: verifier... 
    309             Ch = vkarmn*vkarmn/(func_m*func_h) 
    310             ztmp1 = LOG(zu) - LOG(z0q) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0q*Linv)   ! func_q 
    311             Ce = vkarmn*vkarmn/(func_m*ztmp1) 
    312             ! Non-Solar heat flux to the ocean: 
    313             ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10 
    314             ztmp2 = T_s*T_s 
    315             ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux 
    316                &       +      emiss_w*(rad_lw - stefan*ztmp2*ztmp2)                         ! Net longwave flux 
    317             !!         => "ztmp1" is the net non-solar surface heat flux ! 
    318             !! Updating the values of the skin temperature T_s and q_s : 
    319             CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) 
    320             q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp)  ! 200 -> just to avoid numerics problem on masked regions if silly values are given 
    321          END IF 
    322  
    323          IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN 
     320 
     321         IF( l_use_cs ) THEN 
     322            !! Cool-skin contribution 
     323 
     324            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 
     325               &                   ztmp1, ztmp0,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0 
     326 
     327            CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst, pdTc )  ! Qnsol -> ztmp1 
     328 
     329            T_s(:,:) = zsst(:,:) + pdTc(:,:) 
     330            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:) 
     331            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
     332 
     333         END IF 
     334 
     335         IF( l_use_wl ) THEN 
     336            !! Warm-layer contribution 
     337            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 
     338               &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2 
     339            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst, pdTw ) 
     340            !! Updating T_s and q_s !!! 
     341            T_s(:,:) = zsst(:,:) + pdTw(:,:) 
     342            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:) 
     343            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
     344         END IF 
     345 
     346 
     347         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 
    324348            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
    325349            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    … …  
    337361      Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 
    338362 
    339    END SUBROUTINE TURB_ECMWF 
     363      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc 
     364      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw 
     365 
     366      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
     367      IF (          l_use_cs      ) DEALLOCATE ( pdTc ) 
     368      IF (          l_use_wl      ) DEALLOCATE ( pdTw ) 
     369 
     370   END SUBROUTINE turb_ecmwf 
    340371 
    341372 
    … …  
    433464 
    434465 
    435  
    436  
    437  
    438466   !!====================================================================== 
    439467END MODULE sbcblk_algo_ecmwf 

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90

    r11284 r11615  
    2929   END INTERFACE gamma_moist 
    3030 
     31   INTERFACE e_sat 
     32      MODULE PROCEDURE e_sat_vctr, e_sat_sclr 
     33   END INTERFACE e_sat 
     34 
     35   INTERFACE L_vap 
     36      MODULE PROCEDURE L_vap_vctr, L_vap_sclr 
     37   END INTERFACE L_vap 
     38 
     39   INTERFACE rho_air 
     40      MODULE PROCEDURE rho_air_vctr, rho_air_sclr 
     41   END INTERFACE rho_air 
     42 
     43   INTERFACE cp_air 
     44      MODULE PROCEDURE cp_air_vctr, cp_air_sclr 
     45   END INTERFACE cp_air 
     46 
     47      INTERFACE alpha_sw 
     48      MODULE PROCEDURE alpha_sw_vctr, alpha_sw_sclr 
     49   END INTERFACE alpha_sw 
     50 
     51 
     52 
    3153   PUBLIC virt_temp 
    3254   PUBLIC rho_air 
    … …  
    3961   PUBLIC q_sat 
    4062   PUBLIC q_air_rh 
     63   !: 
     64   PUBLIC update_qnsol_tau 
     65   PUBLIC alpha_sw 
    4166 
    4267   !!---------------------------------------------------------------------- 
    … …  
    7297   END FUNCTION virt_temp 
    7398 
    74    FUNCTION rho_air( ptak, pqa, pslp ) 
    75       !!------------------------------------------------------------------------------- 
    76       !!                           ***  FUNCTION rho_air  *** 
     99   FUNCTION rho_air_vctr( ptak, pqa, pslp ) 
     100      !!------------------------------------------------------------------------------- 
     101      !!                           ***  FUNCTION rho_air_vctr  *** 
    77102      !! 
    78103      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
    … …  
    83108      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
    84109      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
    85       REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
    86       !!------------------------------------------------------------------------------- 
    87       ! 
    88       rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
    89       ! 
    90    END FUNCTION rho_air 
     110      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air_vctr   ! density of moist air   [kg/m^3] 
     111      !!------------------------------------------------------------------------------- 
     112      rho_air_vctr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp ) 
     113   END FUNCTION rho_air_vctr 
     114 
     115   FUNCTION rho_air_sclr( ptak, pqa, pslp ) 
     116      !!------------------------------------------------------------------------------- 
     117      !!                           ***  FUNCTION rho_air_sclr  *** 
     118      !! 
     119      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
     120      !! 
     121      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     122      !!------------------------------------------------------------------------------- 
     123      REAL(wp), INTENT(in) :: ptak           ! air temperature             [K] 
     124      REAL(wp), INTENT(in) :: pqa            ! air specific humidity   [kg/kg] 
     125      REAL(wp), INTENT(in) :: pslp           ! pressure in                [Pa] 
     126      REAL(wp)             :: rho_air_sclr   ! density of moist air   [kg/m^3] 
     127      !!------------------------------------------------------------------------------- 
     128      rho_air_sclr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp ) 
     129   END FUNCTION rho_air_sclr 
     130    
     131 
    91132 
    92133   FUNCTION visc_air(ptak) 
    … …  
    113154   END FUNCTION visc_air 
    114155 
    115    FUNCTION L_vap( psst ) 
     156   FUNCTION L_vap_vctr( psst ) 
    116157      !!--------------------------------------------------------------------------------- 
    117       !!                           ***  FUNCTION L_vap  *** 
    118       !! 
    119       !! ** Purpose : Compute the latent heat of vaporization of water out of temperature 
    120       !! 
    121       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    122       !!---------------------------------------------------------------------------------- 
    123       REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
     158      !!                           ***  FUNCTION L_vap_vctr  *** 
     159      !! 
     160      !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
     161      !! 
     162      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     163      !!---------------------------------------------------------------------------------- 
     164      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap_vctr   ! latent heat of vaporization   [J/kg] 
    124165      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
    125166      !!---------------------------------------------------------------------------------- 
    126167      ! 
    127       L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
    128       ! 
    129    END FUNCTION L_vap 
    130  
    131    FUNCTION cp_air( pqa ) 
    132       !!------------------------------------------------------------------------------- 
    133       !!                           ***  FUNCTION cp_air  *** 
     168      L_vap_vctr = (  2.501_wp - 0.00237_wp * ( psst(:,:) - rt0)  ) * 1.e6_wp 
     169      ! 
     170   END FUNCTION L_vap_vctr 
     171 
     172   FUNCTION L_vap_sclr( psst ) 
     173      !!--------------------------------------------------------------------------------- 
     174      !!                           ***  FUNCTION L_vap_sclr  *** 
     175      !! 
     176      !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
     177      !! 
     178      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     179      !!---------------------------------------------------------------------------------- 
     180      REAL(wp)             ::   L_vap_sclr   ! latent heat of vaporization   [J/kg] 
     181      REAL(wp), INTENT(in) ::   psst         ! water temperature                [K] 
     182      !!---------------------------------------------------------------------------------- 
     183      ! 
     184      L_vap_sclr = (  2.501_wp - 0.00237_wp * ( psst - rt0)  ) * 1.e6_wp 
     185      ! 
     186   END FUNCTION L_vap_sclr 
     187 
     188 
     189   FUNCTION cp_air_vctr( pqa ) 
     190      !!------------------------------------------------------------------------------- 
     191      !!                           ***  FUNCTION cp_air_vctr  *** 
    134192      !! 
    135193      !! ** Purpose : Compute specific heat (Cp) of moist air 
    … …  
    138196      !!------------------------------------------------------------------------------- 
    139197      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
    140       REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
    141       !!------------------------------------------------------------------------------- 
    142       ! 
    143       cp_air = rCp_dry + rCp_vap * pqa 
    144       ! 
    145    END FUNCTION cp_air 
     198      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air_vctr   ! specific heat of moist air   [J/K/kg] 
     199      !!------------------------------------------------------------------------------- 
     200      cp_air_vctr = rCp_dry + rCp_vap * pqa 
     201   END FUNCTION cp_air_vctr 
     202 
     203   FUNCTION cp_air_sclr( pqa ) 
     204      !!------------------------------------------------------------------------------- 
     205      !!                           ***  FUNCTION cp_air_sclr  *** 
     206      !! 
     207      !! ** Purpose : Compute specific heat (Cp) of moist air 
     208      !! 
     209      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     210      !!------------------------------------------------------------------------------- 
     211      REAL(wp), INTENT(in) :: pqa           ! air specific humidity         [kg/kg] 
     212      REAL(wp)             :: cp_air_sclr   ! specific heat of moist air   [J/K/kg] 
     213      !!------------------------------------------------------------------------------- 
     214      cp_air_sclr = rCp_dry + rCp_vap * pqa 
     215   END FUNCTION cp_air_sclr 
     216 
     217 
     218 
     219 
    146220 
    147221   FUNCTION gamma_moist_vctr( ptak, pqa ) 
    … …  
    275349 
    276350 
    277    FUNCTION e_sat_sclr( ptak, pslp ) 
     351   FUNCTION e_sat_vctr(ptak) 
     352      !!************************************************** 
     353      !! ptak:     air temperature [K] 
     354      !! e_sat:  water vapor at saturation [Pa] 
     355      !! 
     356      !! Recommended by WMO 
     357      !! 
     358      !! Goff, J. A., 1957: Saturation pressure of water on the new kelvin 
     359      !! temperature scale. Transactions of the American society of heating 
     360      !! and ventilating engineers, 347–354. 
     361      !! 
     362      !! rt0 should be 273.16 (triple point of water) and not 273.15 like here 
     363      !!************************************************** 
     364 
     365      REAL(wp), DIMENSION(jpi,jpj)             :: e_sat_vctr !: vapour pressure at saturation  [Pa] 
     366      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak    !: temperature (K) 
     367 
     368      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmp 
     369 
     370      ALLOCATE ( ztmp(jpi,jpj) ) 
     371 
     372      ztmp(:,:) = rtt0/ptak(:,:) 
     373 
     374      e_sat_vctr = 100.*( 10.**(10.79574*(1. - ztmp) - 5.028*LOG10(ptak/rtt0)         & 
     375         &       + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak/rtt0 - 1.)) )   & 
     376         &       + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) ) 
     377 
     378      DEALLOCATE ( ztmp ) 
     379 
     380   END FUNCTION e_sat_vctr 
     381 
     382 
     383   FUNCTION e_sat_sclr( ptak ) 
    278384      !!---------------------------------------------------------------------------------- 
    279385      !!                   ***  FUNCTION e_sat_sclr  *** 
    … …  
    287393      !!---------------------------------------------------------------------------------- 
    288394      REAL(wp), INTENT(in) ::   ptak    ! air temperature                  [K] 
    289       REAL(wp), INTENT(in) ::   pslp    ! sea level atmospheric pressure   [Pa] 
    290395      REAL(wp)             ::   e_sat_sclr   ! water vapor at saturation   [kg/kg] 
    291396      ! 
    … …  
    325430         DO ji = 1, jpi 
    326431            ! 
    327             ze_sat =  e_sat_sclr( ptak(ji,jj), pslp(ji,jj) ) 
     432            ze_sat =  e_sat_sclr( ptak(ji,jj) ) 
    328433            ! 
    329434            q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat ) 
    … …  
    351456      DO jj = 1, jpj 
    352457         DO ji = 1, jpi 
    353             ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj), pslp(ji,jj)) 
     458            ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj)) 
    354459            q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze) 
    355460         END DO 
    … …  
    357462      ! 
    358463   END FUNCTION q_air_rh 
     464 
     465 
     466   SUBROUTINE UPDATE_QNSOL_TAU( pTs, pqs, pTa, pqa, pust, ptst, pqst, pUb, pslp, prlw, & 
     467      &                         pQns, pTau,  & 
     468      &                         Qlat) 
     469      !!---------------------------------------------------------------------------------- 
     470      !! Purpose: returns the non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" 
     471      !!          and the module of the wind stress => pTau = Tau 
     472      !! ** Author: L. Brodeau, Sept. 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     473      !!---------------------------------------------------------------------------------- 
     474      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K] 
     475      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg] 
     476      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! air temperature at z=zu [K] 
     477      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=zu [kg/kg] 
     478      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pust ! u* 
     479      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ptst ! t* 
     480      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqst ! q* 
     481      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=zu [m/s] 
     482      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa] 
     483      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: prlw ! downwelling longwave radiative flux [W/m^2] 
     484      ! 
     485      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]] 
     486      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2] 
     487      ! 
     488      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(out) :: Qlat 
     489      ! 
     490      REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zUrho, zTs2, zz0, & 
     491         &        zQlat, zQsen, zQlw 
     492      INTEGER  ::   ji, jj     ! dummy loop indices 
     493      !!----------------------------------------------------------------------------------      
     494      DO jj = 1, jpj 
     495         DO ji = 1, jpi 
     496 
     497            zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt ) 
     498            zdq = pqa(ji,jj) - pqs(ji,jj) ;  zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq ) 
     499            zz0 = pust(ji,jj)/pUb(ji,jj) 
     500            zCd = zz0*zz0 
     501            zCh = zz0*ptst(ji,jj)/zdt 
     502            zCe = zz0*pqst(ji,jj)/zdq 
     503 
     504            zUrho = pUb(ji,jj)*MAX(rho_air(pTa(ji,jj), pqa(ji,jj), pslp(ji,jj)), 1._wp)     ! rho*U10 
     505            zTs2  = pTs(ji,jj)*pTs(ji,jj) 
     506 
     507            ! Wind stress module: 
     508            pTau(ji,jj) = zCd*zUrho*pUb(ji,jj) ! lolo? 
     509 
     510            ! Non-Solar heat flux to the ocean: 
     511            zQlat = MIN ( zUrho*zCe*L_vap( pTs(ji,jj)) * zdq , 0._wp )  ! we do not want a Qlat > 0 ! 
     512            zQsen = zUrho*zCh*cp_air(pqa(ji,jj)) * zdt 
     513            zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux 
     514 
     515            pQns(ji,jj) = zQlat + zQsen + zQlw 
     516             
     517            IF ( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat             
     518         END DO 
     519      END DO 
     520   END SUBROUTINE UPDATE_QNSOL_TAU 
     521 
     522 
     523   FUNCTION alpha_sw_vctr( psst ) 
     524      !!--------------------------------------------------------------------------------- 
     525      !!                           ***  FUNCTION alpha_sw_vctr  *** 
     526      !! 
     527      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa) 
     528      !! 
     529      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     530      !!---------------------------------------------------------------------------------- 
     531      REAL(wp), DIMENSION(jpi,jpj)             ::   alpha_sw_vctr   ! latent heat of vaporization   [J/kg] 
     532      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
     533      !!---------------------------------------------------------------------------------- 
     534      alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79 
     535   END FUNCTION alpha_sw_vctr 
     536 
     537   FUNCTION alpha_sw_sclr( psst ) 
     538      !!--------------------------------------------------------------------------------- 
     539      !!                           ***  FUNCTION alpha_sw_sclr  *** 
     540      !! 
     541      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa) 
     542      !! 
     543      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     544      !!---------------------------------------------------------------------------------- 
     545      REAL(wp)             ::   alpha_sw_sclr   ! latent heat of vaporization   [J/kg] 
     546      REAL(wp), INTENT(in) ::   psst   ! sea-water temperature                   [K] 
     547      !!---------------------------------------------------------------------------------- 
     548      alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79 
     549   END FUNCTION alpha_sw_sclr 
     550 
    359551 
    360552 

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/demo_cfgs.txt

    r10516 r11615  
    99WAD OCE 
    1010BENCH OCE ICE TOP 
     11STATION_TASF OCE 
Note: See TracChangeset for help on using the changeset viewer.