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 11772 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src – NEMO

Ignore:
Timestamp:
2019-10-23T16:04:12+02:00 (4 years ago)
Author:
laurent
Message:

LB: solid updates+improvements of cool-skin/warm-layer capabilty of COARE and ECMWF bulk algorithms!

Location:
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE
Files:
8 edited

Legend:

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

    r11615 r11772  
    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 
     84   REAL(wp), PARAMETER, PUBLIC :: radrw = rho0_a/rho0_w !: Density ratio 
     85   REAL(wp), PARAMETER, PUBLIC :: sq_radrw = SQRT(rho0_a/rho0_w) 
    8586   REAL(wp), PARAMETER, PUBLIC :: rCp0_w  = 4190._wp    !: Specific heat capacity of seawater (ECMWF 4190) [J/K/kg] 
    8687   REAL(wp), PARAMETER, PUBLIC :: rnu0_w  = 1.e-6_wp    !: kinetic viscosity of water                      [m^2/s] 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90

    r11672 r11772  
    542542         CASE( np_ECMWF     ) 
    543543            IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_ecmwf" WITH CSWL options!!!, gdept_1d(1)=', gdept_1d(1) !LBrm 
    544             CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF 
     544            CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF 
    545545               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
    546546               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
     
    585585         CASE( np_ECMWF     ) 
    586586            IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!' !LBrm 
    587             CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF 
     587            CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF 
    588588               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    589589 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90

    r11672 r11772  
    4949   REAL(wp), PARAMETER ::   Beta0   =   1.25_wp   ! gustiness parameter 
    5050 
    51    INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations 
     51   INTEGER , PARAMETER ::   nb_itt = 10             ! number of itterations 
    5252 
    5353   !!---------------------------------------------------------------------- 
     
    7070      IF ( l_use_wl ) THEN 
    7171         ierr = 0 
    72          ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), H_wl(jpi,jpj), STAT=ierr ) 
     72         ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), Hz_wl(jpi,jpj), dT_wl(jpi,jpj), STAT=ierr ) 
    7373         !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of Tau_ac and Qnt_ac failed!' 
    7474         Tau_ac(:,:) = 0._wp 
    75          Qnt_ac(:,:)   = 0._wp 
    76          H_wl(:,:)    = H_wl_max 
     75         Qnt_ac(:,:) = 0._wp 
     76         Hz_wl(:,:)  = Hwl_max 
     77         dT_wl(:,:)  = 0._wp 
    7778      END IF 
    7879      !! 
    7980      IF ( l_use_cs ) THEN 
    8081         ierr = 0 
    81          ALLOCATE ( delta_vl(jpi,jpj), STAT=ierr ) 
    82          !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of delta_vl and Qnt_ac failed!' 
    83          delta_vl(:,:) = 0.001_wp      ! First guess of zdelta [m] 
     82         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 
     83         !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of dT_cs and Qnt_ac failed!' 
     84         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction 
    8485      END IF 
    8586   END SUBROUTINE coare3p0_init 
     
    9192      &                      Cdn, Chn, Cen,                                              & 
    9293      &                      Qsw, rad_lw, slp, pdT_cs,                                   & ! optionals for cool-skin (and warm-layer) 
    93       &                      pdT_wl, Hwl )                                                 ! optionals for warm-layer only 
     94      &                      pdT_wl, pHz_wl )                                                 ! optionals for warm-layer only 
    9495      !!---------------------------------------------------------------------- 
    9596      !!                      ***  ROUTINE  turb_coare3p0  *** 
     
    133134      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K] 
    134135      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K] 
    135       !!    * Hwl     : depth of warm layer                                   [m] 
     136      !!    * pHz_wl  : thickness of warm-layer                               [m] 
    136137      !! 
    137138      !! OUTPUT : 
     
    171172      ! 
    172173      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K] 
    173       REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Hwl      !             [m] 
     174      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m] 
    174175      ! 
    175176      INTEGER :: j_itt 
     
    190191         &                pdTw,   &  ! SST increment "dT" for warm layer correction          [K] 
    191192         &                zHwl       ! depth of warm-layer [m] 
     193      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0' 
    192194      !!---------------------------------------------------------------------------------- 
    193195 
    194       IF ( kt == 1 ) CALL COARE3P0_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays 
     196      IF ( kt == nit000 ) CALL COARE3P0_INIT(l_use_cs, l_use_wl) 
    195197 
    196198      l_zt_equal_zu = .FALSE. 
     
    199201 
    200202      !! Initializations for cool skin and warm layer: 
     203      IF ( l_use_cs ) THEN 
     204         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 
     205            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'; STOP 
     206         END IF 
     207      END IF 
     208 
     209      IF ( l_use_wl ) THEN 
     210         IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) )) THEN 
     211            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!'; STOP 
     212         END IF 
     213      END IF 
     214 
    201215      IF ( l_use_cs .OR. l_use_wl ) THEN 
    202          IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) & 
    203             & CALL ctl_stop( 'turb_coare3p0 => provide Qsw, rad_lw & slp to ', 'use cool-skin and/or warm-layer param' ) 
    204216         ALLOCATE ( zsst(jpi,jpj) ) 
    205217         zsst = T_s ! backing up the bulk SST 
    206          IF( l_use_cs ) T_s = T_s - 0.25   ! First guess of correction 
     218         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction 
    207219         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!! 
    208       END IF 
    209       IF ( l_use_cs ) THEN 
    210          ALLOCATE ( pdTc(jpi,jpj) ) 
    211          pdTc(:,:) = -0.25_wp  ! First guess of skin correction 
    212       END IF 
    213       IF ( l_use_wl ) THEN 
    214          ALLOCATE ( pdTw(jpi,jpj) ) 
    215          IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) ) 
    216220      END IF 
    217221 
     
    323327            !! Cool-skin contribution 
    324328 
    325             CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 
     329            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 
    326330               &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
    327331 
    328             CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2,  pdTc )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2 
    329  
    330             T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1) 
    331             IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)*tmask(:,:,1) 
     332            CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2 
     333 
     334            T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1) 
     335            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 
    332336            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    333337 
     
    336340         IF( l_use_wl ) THEN 
    337341            !! Warm-layer contribution 
    338             CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 
     342            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 
    339343               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
    340344            !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 
    341             IF (PRESENT(Hwl)) THEN 
    342                CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl ) 
    343             ELSE 
    344                CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt),  pdTw ) 
    345             END IF 
     345            CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt) ) 
     346 
    346347            !! Updating T_s and q_s !!! 
    347             T_s(:,:) = zsst(:,:) + pdTw(:,:)*tmask(:,:,1) 
    348             IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)*tmask(:,:,1) 
     348            T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) 
     349            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 
    349350            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    350             !LOLO: 
    351             !IF ( j_itt == nb_itt) THEN 
    352             !   WRITE(cf_tmp,'("Qnt_ac_k",i5.5,"_p",i4.4,".nc")') kt, narea 
    353             !   CALL DUMP_FIELD(REAL( Qnt_ac*tmask(:,:,1) , 4), TRIM(cf_tmp), 'Qnt_ac') 
    354             !   WRITE(cf_tmp,  '("pdTw_k",i5.5,"_p",i4.4,".nc")') kt, narea 
    355             !   CALL DUMP_FIELD(REAL( pdTw*tmask(:,:,1) , 4), TRIM(cf_tmp), 'pdTw') 
    356             !END IF 
    357             !LOLO. 
    358          END IF 
    359  
     351         END IF 
    360352 
    361353         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 
     
    379371      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 
    380372 
    381       IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc*tmask(:,:,1) 
    382       IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw*tmask(:,:,1) 
     373      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
     374      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 
     375      IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 
    383376 
    384377      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
    385       IF (          l_use_cs      ) DEALLOCATE ( pdTc ) 
    386       IF (          l_use_wl      ) THEN 
    387          DEALLOCATE ( pdTw ) 
    388          IF (PRESENT(Hwl)) DEALLOCATE ( zHwl ) 
    389       END IF 
    390378 
    391379   END SUBROUTINE turb_coare3p0 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90

    r11672 r11772  
    4545   PUBLIC :: COARE3P6_INIT, TURB_COARE3P6 
    4646 
    47    !                                              !! COARE own values for given constants: 
    48    REAL(wp), PARAMETER ::   zi0     = 600._wp     ! scale height of the atmospheric boundary layer... 
    49    REAL(wp), PARAMETER ::   Beta0   =   1.2_wp   ! gustiness parameter 
    50  
    51    INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations 
     47   !! COARE own values for given constants: 
     48   REAL(wp), PARAMETER :: zi0   = 600._wp     ! scale height of the atmospheric boundary layer... 
     49   REAL(wp), PARAMETER :: Beta0 =  1.2_wp     ! gustiness parameter 
     50 
     51   INTEGER , PARAMETER ::   nb_itt = 10             ! number of itterations 
    5252 
    5353   !!---------------------------------------------------------------------- 
     
    7070      IF ( l_use_wl ) THEN 
    7171         ierr = 0 
    72          ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), H_wl(jpi,jpj), STAT=ierr ) 
    73          !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of Tau_ac and Qnt_ac failed!' 
     72         ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), Hz_wl(jpi,jpj), dT_wl(jpi,jpj), STAT=ierr ) 
     73         !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of Tau_ac, Qnt_ac, dT_wl & Hz_wl failed!' 
    7474         Tau_ac(:,:) = 0._wp 
    75          Qnt_ac(:,:)   = 0._wp 
    76          H_wl(:,:)    = H_wl_max 
     75         Qnt_ac(:,:) = 0._wp 
     76         dT_wl(:,:)  = 0._wp 
     77         Hz_wl(:,:)  = Hwl_max 
    7778      END IF 
    7879      !! 
    7980      IF ( l_use_cs ) THEN 
    8081         ierr = 0 
    81          ALLOCATE ( delta_vl(jpi,jpj), STAT=ierr ) 
    82          !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of delta_vl and Qnt_ac failed!' 
    83          delta_vl(:,:) = 0.001_wp      ! First guess of zdelta [m] 
     82         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 
     83         !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of dT_cs failed!' 
     84         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction 
    8485      END IF 
    8586   END SUBROUTINE coare3p6_init 
     
    9192      &                      Cdn, Chn, Cen,                                              & 
    9293      &                      Qsw, rad_lw, slp, pdT_cs,                                   & ! optionals for cool-skin (and warm-layer) 
    93       &                      pdT_wl, Hwl )                                                 ! optionals for warm-layer only 
     94      &                      pdT_wl, pHz_wl )                                                 ! optionals for warm-layer only 
    9495      !!---------------------------------------------------------------------- 
    9596      !!                      ***  ROUTINE  turb_coare3p6  *** 
     
    131132      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2] 
    132133      !!    *  slp    : sea-level pressure                                    [Pa] 
     134      !! 
     135      !! OPTIONAL OUTPUT: 
     136      !! ---------------- 
    133137      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K] 
    134138      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K] 
    135       !!    * Hwl     : depth of warm layer                                   [m] 
     139      !!    * pHz_wl  : thickness of warm-layer                               [m] 
    136140      !! 
    137141      !! OUTPUT : 
     
    169173      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa] 
    170174      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs 
    171       ! 
    172175      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K] 
    173       REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Hwl      !             [m] 
     176      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m] 
    174177      ! 
    175178      INTEGER :: j_itt 
     
    190193         &                pdTw,   &  ! SST increment "dT" for warm layer correction          [K] 
    191194         &                zHwl       ! depth of warm-layer [m] 
     195      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p6@sbcblk_algo_coare3p6' 
    192196      !!---------------------------------------------------------------------------------- 
    193197 
    194       IF ( kt == 1 ) CALL COARE3P6_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays 
     198      IF ( kt == nit000 ) CALL COARE3P6_INIT(l_use_cs, l_use_wl) 
    195199 
    196200      l_zt_equal_zu = .FALSE. 
     
    199203 
    200204      !! Initializations for cool skin and warm layer: 
     205      IF ( l_use_cs ) THEN 
     206         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 
     207            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'; STOP 
     208         END IF 
     209      END IF 
     210 
     211      IF ( l_use_wl ) THEN 
     212         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 
     213            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!'; STOP 
     214         END IF 
     215      END IF 
     216 
    201217      IF ( l_use_cs .OR. l_use_wl ) THEN 
    202          IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) & 
    203             & CALL ctl_stop( 'turb_coare3p6 => provide Qsw, rad_lw & slp to ', 'use cool-skin and/or warm-layer param' ) 
    204218         ALLOCATE ( zsst(jpi,jpj) ) 
    205219         zsst = T_s ! backing up the bulk SST 
    206          IF( l_use_cs ) T_s = T_s - 0.25   ! First guess of correction 
     220         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction 
    207221         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!! 
    208       END IF 
    209       IF ( l_use_cs ) THEN 
    210          ALLOCATE ( pdTc(jpi,jpj) ) 
    211          pdTc(:,:) = -0.25_wp  ! First guess of skin correction 
    212       END IF 
    213       IF ( l_use_wl ) THEN 
    214          ALLOCATE ( pdTw(jpi,jpj) ) 
    215          IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) ) 
    216222      END IF 
    217223 
     
    323329            !! Cool-skin contribution 
    324330 
    325             CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 
     331            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 
    326332               &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
    327333 
    328             CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2,  pdTc )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2 
    329  
    330             T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1) 
    331             IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)*tmask(:,:,1) 
     334            CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2 
     335 
     336            T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1) 
     337            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 
    332338            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    333339 
     
    336342         IF( l_use_wl ) THEN 
    337343            !! Warm-layer contribution 
    338             CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 
     344            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 
    339345               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
    340346            !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 
    341             IF (PRESENT(Hwl)) THEN 
    342                CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl ) 
    343             ELSE 
    344                CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt),  pdTw ) 
    345             END IF 
     347            CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt) ) 
     348 
    346349            !! Updating T_s and q_s !!! 
    347             T_s(:,:) = zsst(:,:) + pdTw(:,:)*tmask(:,:,1) 
    348             IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)*tmask(:,:,1) 
     350            T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) 
     351            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 
    349352            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    350             !LOLO: 
    351             !IF ( j_itt == nb_itt) THEN 
    352             !   WRITE(cf_tmp,'("Qnt_ac_k",i5.5,"_p",i4.4,".nc")') kt, narea 
    353             !   CALL DUMP_FIELD(REAL( Qnt_ac*tmask(:,:,1) , 4), TRIM(cf_tmp), 'Qnt_ac') 
    354             !   WRITE(cf_tmp,  '("pdTw_k",i5.5,"_p",i4.4,".nc")') kt, narea 
    355             !   CALL DUMP_FIELD(REAL( pdTw*tmask(:,:,1) , 4), TRIM(cf_tmp), 'pdTw') 
    356             !END IF 
    357             !LOLO. 
    358          END IF 
    359  
     353         END IF 
    360354 
    361355         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 
     
    379373      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 
    380374 
    381       IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc*tmask(:,:,1) 
    382       IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw*tmask(:,:,1) 
     375      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
     376      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 
     377      IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 
    383378 
    384379      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
    385       IF (          l_use_cs      ) DEALLOCATE ( pdTc ) 
    386       IF (          l_use_wl      ) THEN 
    387          DEALLOCATE ( pdTw ) 
    388          IF (PRESENT(Hwl)) DEALLOCATE ( zHwl ) 
    389       END IF 
    390380 
    391381   END SUBROUTINE turb_coare3p6 
     
    395385      !!------------------------------------------------------------------- 
    396386      !! Computes the Charnock parameter as a function of the Neutral wind speed at 10m 
    397       !! 
     387      !!  "wind speed dependent formulation" 
    398388      !!  (Eq. 13 in Edson et al., 2013) 
    399389      !! 
     
    405395      REAL(wp), PARAMETER :: charn0_max = 0.028  !: value above which the Charnock parameter levels off for winds > 18 m/s 
    406396      !!------------------------------------------------------------------- 
    407       ! 
    408397      alfa_charn_3p6 = MAX( MIN( 0.0017_wp*pwnd - 0.005_wp , charn0_max) , 0._wp ) 
    409       ! 
     398      !! 
    410399   END FUNCTION alfa_charn_3p6 
     400 
     401   FUNCTION alfa_charn_3p6_wave( pus, pwsh, pwps ) 
     402      !!------------------------------------------------------------------- 
     403      !! Computes the Charnock parameter as a function of wave information and u* 
     404      !! 
     405      !!  (COARE 3.6, Fairall et al., 2018) 
     406      !! 
     407      !! Author: L. Brodeau, October 2019 / AeroBulk  (https://github.com/brodeau/aerobulk/) 
     408      !!------------------------------------------------------------------- 
     409      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p6_wave 
     410      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus   ! friction velocity             [m/s] 
     411      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwsh  ! significant wave height       [m] 
     412      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwps  ! phase speed of dominant waves [m/s] 
     413      !!------------------------------------------------------------------- 
     414      alfa_charn_3p6_wave = ( pwsh*0.2_wp*(pus/pwps)**2.2_wp ) * grav/(pus*pus) 
     415      !! 
     416   END FUNCTION alfa_charn_3p6_wave 
     417 
    411418 
    412419   FUNCTION psi_m_coare( pzeta ) 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r11631 r11772  
    4444   PRIVATE 
    4545 
    46    PUBLIC :: TURB_ECMWF   ! called by sbcblk.F90 
     46   PUBLIC :: ECMWF_INIT, TURB_ECMWF 
    4747 
    4848   !                   !! ECMWF own values for given constants, taken form IFS documentation... 
     
    5555   REAL(wp), PARAMETER ::   alpha_Q = 0.62    ! 
    5656 
    57    INTEGER , PARAMETER ::   nb_itt = 5        ! number of itterations 
     57   INTEGER , PARAMETER ::   nb_itt = 10             ! number of itterations 
    5858 
    5959   !!---------------------------------------------------------------------- 
    6060CONTAINS 
    6161 
    62    SUBROUTINE turb_ecmwf( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 
    63       &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
    64       &                      Cdn, Chn, Cen,                      & 
    65       &                      Qsw, rad_lw, slp, pdT_cs,           & ! optionals for cool-skin (and warm-layer) 
    66       &                      pdT_wl )                              ! optionals for warm-layer only 
     62 
     63   SUBROUTINE ecmwf_init(l_use_cs, l_use_wl) 
     64      !!--------------------------------------------------------------------- 
     65      !!                  ***  FUNCTION ecmwf_init  *** 
     66      !! 
     67      !! INPUT : 
     68      !! ------- 
     69      !!    * l_use_cs : use the cool-skin parameterization 
     70      !!    * l_use_wl : use the warm-layer parameterization 
     71      !!--------------------------------------------------------------------- 
     72      LOGICAL , INTENT(in) ::   l_use_cs ! use the cool-skin parameterization 
     73      LOGICAL , INTENT(in) ::   l_use_wl ! use the warm-layer parameterization 
     74      INTEGER :: ierr 
     75      !!--------------------------------------------------------------------- 
     76      IF ( l_use_wl ) THEN 
     77         ierr = 0 
     78         ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) 
     79         !IF( ierr > 0 ) STOP ' ECMWF_INIT => allocation of dT_wl & Hz_wl failed!' 
     80         dT_wl(:,:)  = 0._wp 
     81         Hz_wl(:,:)  = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars) 
     82      END IF 
     83      !! 
     84      IF ( l_use_cs ) THEN 
     85         ierr = 0 
     86         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 
     87         !IF( ierr > 0 ) STOP ' ECMWF_INIT => allocation of dT_cs failed!' 
     88         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction 
     89      END IF 
     90   END SUBROUTINE ecmwf_init 
     91 
     92 
     93 
     94   SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 
     95      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                           & 
     96      &                      Cdn, Chn, Cen,                                           & 
     97      &                      Qsw, rad_lw, slp, pdT_cs,                                & ! optionals for cool-skin (and warm-layer) 
     98      &                      pdT_wl, pHz_wl )                                                 ! optionals for warm-layer only 
    6799      !!---------------------------------------------------------------------- 
    68100      !!                      ***  ROUTINE  turb_ecmwf  *** 
     
    80112      !! INPUT : 
    81113      !! ------- 
     114      !!    *  kt   : current time step (starts at 1) 
    82115      !!    *  zt   : height for temperature and spec. hum. of air            [m] 
    83116      !!    *  zu   : height for wind speed (usually 10m)                     [m] 
     
    95128      !! 
    96129      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg] 
    97       !!              -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True) 
    98       !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False) 
     130      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True) 
     131      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False) 
    99132      !! 
    100133      !! OPTIONAL INPUT: 
     
    103136      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2] 
    104137      !!    *  slp    : sea-level pressure                                    [Pa] 
     138      !! 
     139      !! OPTIONAL OUTPUT: 
     140      !! ---------------- 
    105141      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K] 
    106142      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K] 
     143      !!    * pHz_wl  : thickness of warm-layer                               [m] 
    107144      !! 
    108145      !! OUTPUT : 
     
    118155      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    119156      !!---------------------------------------------------------------------------------- 
     157      INTEGER,  INTENT(in   )                     ::   kt       ! current time step 
    120158      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
    121159      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m] 
     
    139177      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa] 
    140178      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs 
    141       ! 
    142179      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K] 
     180      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m] 
    143181      ! 
    144182      INTEGER :: j_itt 
     
    162200      !!---------------------------------------------------------------------------------- 
    163201 
     202      IF ( kt == nit000 ) CALL ECMWF_INIT(l_use_cs, l_use_wl) 
     203 
    164204      l_zt_equal_zu = .FALSE. 
    165205      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     
    168208      IF ( l_use_cs ) THEN 
    169209         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 
     210            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'; STOP 
     211         END IF 
    175212      END IF 
    176213 
    177214      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 
     215         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) THEN 
     216            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!'; STOP 
     217         END IF 
    184218      END IF 
    185219 
     
    322356            !! Cool-skin contribution 
    323357 
    324             CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 
     358            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 
    325359               &                   ztmp1, ztmp0,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0 
    326360 
    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(:,:) 
     361            CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst )  ! Qnsol -> ztmp1 
     362 
     363            T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1) 
     364            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 
    331365            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    332366 
     
    335369         IF( l_use_wl ) THEN 
    336370            !! 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, & 
     371            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 
    338372               &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2 
    339             CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst, pdTw ) 
     373            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst ) 
    340374            !! Updating T_s and q_s !!! 
    341             T_s(:,:) = zsst(:,:) + pdTw(:,:) 
    342             IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:) 
     375            T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) 
     376            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 
    343377            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    344378         END IF 
    345  
    346379 
    347380         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 
     
    361394      Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 
    362395 
    363       IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc 
    364       IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw 
     396      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
     397      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 
     398      IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 
    365399 
    366400      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
    367       IF (          l_use_cs      ) DEALLOCATE ( pdTc ) 
    368       IF (          l_use_wl      ) DEALLOCATE ( pdTw ) 
    369401 
    370402   END SUBROUTINE turb_ecmwf 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90

    r11638 r11772  
    4545   END INTERFACE cp_air 
    4646 
    47       INTERFACE alpha_sw 
     47   INTERFACE alpha_sw 
    4848      MODULE PROCEDURE alpha_sw_vctr, alpha_sw_sclr 
    4949   END INTERFACE alpha_sw 
     50 
     51   INTERFACE turb_fluxes 
     52      MODULE PROCEDURE turb_fluxes_vctr, turb_fluxes_sclr 
     53   END INTERFACE turb_fluxes 
    5054 
    5155 
     
    6468   PUBLIC update_qnsol_tau 
    6569   PUBLIC alpha_sw 
     70   PUBLIC turb_fluxes 
    6671 
    6772   !!---------------------------------------------------------------------- 
     
    128133      rho_air_sclr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp ) 
    129134   END FUNCTION rho_air_sclr 
    130     
     135 
    131136 
    132137 
     
    464469 
    465470 
    466    SUBROUTINE UPDATE_QNSOL_TAU( pTs, pqs, pTa, pqa, pust, ptst, pqst, pUb, pslp, prlw, & 
     471   SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, pslp, prlw, & 
    467472      &                         pQns, pTau,  & 
    468473      &                         Qlat) 
     
    472477      !! ** Author: L. Brodeau, Sept. 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    473478      !!---------------------------------------------------------------------------------- 
     479      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m) 
    474480      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K] 
    475481      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] 
     482      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K] 
     483      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg] 
    478484      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pust ! u* 
    479485      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ptst ! t* 
    480486      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] 
     487      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s] 
     488      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 
    482489      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa] 
    483490      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: prlw ! downwelling longwave radiative flux [W/m^2] 
     
    491498         &        zQlat, zQsen, zQlw 
    492499      INTEGER  ::   ji, jj     ! dummy loop indices 
    493       !!----------------------------------------------------------------------------------      
     500      !!---------------------------------------------------------------------------------- 
    494501      DO jj = 1, jpj 
    495502         DO ji = 1, jpi 
     
    502509            zCe = zz0*pqst(ji,jj)/zdq 
    503510 
    504             zUrho = pUb(ji,jj)*MAX(rho_air(pTa(ji,jj), pqa(ji,jj), pslp(ji,jj)), 1._wp)     ! rho*U10 
     511            !zUrho = pUb(ji,jj)*MAX(rho_air(pTa(ji,jj), pqa(ji,jj), pslp(ji,jj)), 1._wp)     ! rho*U10 
    505512            zTs2  = pTs(ji,jj)*pTs(ji,jj) 
    506513 
     514            CALL TURB_FLUXES( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, & 
     515               &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj), & 
     516               &              pTau(ji,jj), zQsen, zQlat ) 
     517 
     518 
    507519            ! Wind stress module: 
    508             pTau(ji,jj) = zCd*zUrho*pUb(ji,jj) ! lolo? 
     520            !pTau(ji,jj) = zCd*zUrho*pUb(ji,jj) ! lolo? 
    509521 
    510522            ! 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 
     523            !zQlat = MIN ( zUrho*zCe*L_vap( pTs(ji,jj)) * zdq , 0._wp )  ! we do not want a Qlat > 0 ! 
     524            !zQsen = zUrho*zCh*cp_air(pqa(ji,jj)) * zdt 
    513525            zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux 
    514526 
    515527            pQns(ji,jj) = zQlat + zQsen + zQlw 
    516              
    517             IF ( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat             
     528 
     529            IF ( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 
    518530         END DO 
    519531      END DO 
    520532   END SUBROUTINE UPDATE_QNSOL_TAU 
     533 
     534 
     535 
     536 
     537 
     538   SUBROUTINE TURB_FLUXES_VCTR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, & 
     539      &                                 pTau, pQsen, pQlat,  pEvap ) 
     540      !!---------------------------------------------------------------------------------- 
     541      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m) 
     542      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K] 
     543      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg] 
     544      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K] 
     545      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg] 
     546      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCd 
     547      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCh 
     548      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCe 
     549      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s] 
     550      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 
     551      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa] 
     552      !! 
     553      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau  ! module of the wind stress [N/m^2] 
     554      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQsen !  [W/m^2] 
     555      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat !  [W/m^2] 
     556      !! 
     557      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap !  [kg/m^2/s] 
     558      !! 
     559      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap 
     560      INTEGER  :: ji, jj, jq     ! dummy loop indices 
     561      !!---------------------------------------------------------------------------------- 
     562      DO jj = 1, jpj 
     563         DO ji = 1, jpi 
     564 
     565            !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 
     566            ztaa = pTa(ji,jj) ! first guess... 
     567            DO jq = 1, 4 
     568               zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) ) 
     569               ztaa = pTa(ji,jj) - zgamma*pzu   ! Absolute temp. is slightly colder... 
     570            END DO 
     571            zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)) 
     572            zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 
     573 
     574            zUrho = pUb(ji,jj)*MAX(zrho, 1._wp)     ! rho*U10 
     575 
     576            pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module 
     577 
     578            zevap        = MIN( zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj)) , 0._wp )   ! we do not want condensation & Qlat > 0 ! 
     579            pQsen(ji,jj) =      zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj)) 
     580            pQlat(ji,jj) =  L_vap(pTs(ji,jj)) * zevap 
     581 
     582            IF ( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap 
     583 
     584         END DO 
     585      END DO 
     586   END SUBROUTINE TURB_FLUXES_VCTR 
     587 
     588 
     589   SUBROUTINE TURB_FLUXES_SCLR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, & 
     590      &                                 pTau, pQsen, pQlat,  pEvap ) 
     591      !!---------------------------------------------------------------------------------- 
     592      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m) 
     593      REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K] 
     594      REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg] 
     595      REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K] 
     596      REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg] 
     597      REAL(wp), INTENT(in)  :: pCd 
     598      REAL(wp), INTENT(in)  :: pCh 
     599      REAL(wp), INTENT(in)  :: pCe 
     600      REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s] 
     601      REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 
     602      REAL(wp), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa] 
     603      !! 
     604      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2] 
     605      REAL(wp), INTENT(out) :: pQsen !  [W/m^2] 
     606      REAL(wp), INTENT(out) :: pQlat !  [W/m^2] 
     607      !! 
     608      REAL(wp), INTENT(out), OPTIONAL :: pEvap !  [kg/m^2/s] 
     609      !! 
     610      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap 
     611      INTEGER  :: jq 
     612      !!---------------------------------------------------------------------------------- 
     613 
     614      !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 
     615      ztaa = pTa ! first guess... 
     616      DO jq = 1, 4 
     617         zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa ) 
     618         ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder... 
     619      END DO 
     620      zrho = rho_air(ztaa, pqa, pslp) 
     621      zrho = rho_air(ztaa, pqa, pslp-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 
     622 
     623      zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10 
     624 
     625      pTau = zUrho * pCd * pwnd ! Wind stress module 
     626 
     627      zevap        = MIN( zUrho * pCe * (pqa - pqs) , 0._wp )   ! we do not want condensation & Qlat > 0 ! 
     628      pQsen =      zUrho * pCh * (pTa - pTs) * cp_air(pqa) 
     629      pQlat =  L_vap(pTs) * zevap 
     630 
     631      IF ( PRESENT(pEvap) ) pEvap = - zevap 
     632 
     633   END SUBROUTINE TURB_FLUXES_SCLR 
     634 
     635 
     636 
    521637 
    522638 
     
    529645      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    530646      !!---------------------------------------------------------------------------------- 
    531       REAL(wp), DIMENSION(jpi,jpj)             ::   alpha_sw_vctr   ! latent heat of vaporization   [J/kg] 
     647      REAL(wp), DIMENSION(jpi,jpj)             ::   alpha_sw_vctr   ! thermal expansion coefficient of sea-water [1/K] 
    532648      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
    533649      !!---------------------------------------------------------------------------------- 
     
    543659      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    544660      !!---------------------------------------------------------------------------------- 
    545       REAL(wp)             ::   alpha_sw_sclr   ! latent heat of vaporization   [J/kg] 
     661      REAL(wp)             ::   alpha_sw_sclr   ! thermal expansion coefficient of sea-water [1/K] 
    546662      REAL(wp), INTENT(in) ::   psst   ! sea-water temperature                   [K] 
    547663      !!---------------------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90

    r11672 r11772  
    3838 
    3939   !! Cool-skin related parameters: 
    40    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: delta_vl  !: thickness of the surface viscous layer (in the water) right below the air-sea interface [m] 
     40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 
     41      &                        dT_cs         !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta) 
     42   REAL(wp), PARAMETER :: zcon0 = -16._wp * grav * rho0_w * rCp0_w * rnu0_w*rnu0_w*rnu0_w / ( rk0_w*rk0_w ) ! "-" because ocean convention: Qabs > 0 => gain of heat for ocean! 
     43   !!                             => see eq.(14) in Fairall et al. 1996   (eq.(6) of Zeng aand Beljaars is WRONG! (typo?) 
    4144 
    4245   !! Warm-layer related parameters: 
    43    REAL(wp), PARAMETER, PUBLIC :: H_wl_max = 20._wp    !: maximum depth of warm layer (adjustable) 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 
     47      &                        dT_wl,     &  !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 
     48      &                        Hz_wl,     &  !: depth of warm-layer [m] 
     49      &                        Qnt_ac,    &  !: time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight) 
     50      &                        Tau_ac        !: time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight) 
     51 
     52   REAL(wp), PARAMETER, PUBLIC :: Hwl_max = 20._wp    !: maximum depth of warm layer (adjustable) 
    4453   ! 
    4554   REAL(wp), PARAMETER :: rich   = 0.65_wp   !: critical Richardson number 
    4655   ! 
    47    REAL(wp), PARAMETER :: Qabs_thr = 50._wp  !: threshold for heat flux absorbed in WL 
    4856   REAL(wp), PARAMETER :: zfr0   = 0.5_wp    !: initial value of solar flux absorption 
    4957   ! 
    50    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Tau_ac  ! time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight) 
    51    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Qnt_ac    ! time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight) 
    52    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: H_wl     ! depth of warm-layer [m] 
    5358   !!---------------------------------------------------------------------- 
    5459CONTAINS 
    5560 
    5661 
    57    SUBROUTINE CS_COARE( pQsw, pQnsol, pustar, pSST, pQlat,  pdT ) 
     62   SUBROUTINE CS_COARE( pQsw, pQnsol, pustar, pSST, pQlat ) 
    5863      !!--------------------------------------------------------------------- 
    5964      !! 
    6065      !!  Cool-Skin scheme according to Fairall et al. 1996, revisited for COARE 3.6 (Fairall et al., 2019) 
    61       !!     ------------------------------------------------------------------ 
     66      !! 
     67      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 
     68      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 
     69      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 
     70      !! doi:10.1029/95JC03190. 
     71      !! 
     72      !!------------------------------------------------------------------ 
    6273      !! 
    6374      !!  **   INPUT: 
     
    6879      !!     *pQlat*      surface latent heat flux                       [K] 
    6980      !! 
    70       !!  **  INPUT/OUTPUT: 
    71       !!     *pdT*  : as input =>  previous estimate of dT cool-skin 
    72       !!              as output =>  new estimate of dT cool-skin 
    73       !! 
    7481      !!------------------------------------------------------------------ 
    75       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 
    76       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 
    77       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*) 
    78       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pSST ! bulk SST [K] 
    79       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQlat  ! latent heat flux [W/m^2] 
    80       !! 
    81       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT    ! dT due to cool-skin effect 
    82       !!--------------------------------------------------------------------- 
    83       INTEGER  ::   ji, jj     ! dummy loop indices 
    84       REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, & 
    85          &        zz1, zz2, zus, & 
    86          &        ztf 
    87       !!--------------------------------------------------------------------- 
    88  
     82      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 
     83      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 
     84      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*) 
     85      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST ! bulk SST [K] 
     86      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQlat  ! latent heat flux [W/m^2] 
     87      !!--------------------------------------------------------------------- 
     88      INTEGER  :: ji, jj, jc 
     89      REAL(wp) :: zQabs, zdelta, zfr 
     90      !!--------------------------------------------------------------------- 
    8991      DO jj = 1, jpj 
    9092         DO ji = 1, jpi 
    9193 
    92             zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 
    93  
    94             zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere 
    95  
    96             zdelta = delta_vl(ji,jj)   ! using last value of delta 
    97  
    98             !! Fraction of the shortwave flux absorbed by the cool-skin sublayer: 
    99             zfr   = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Eq.16 (Fairall al. 1996b) /  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 ! 
    100  
    101             zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) )  ! Total cooling at the interface 
    102  
    103             ztf = 0.5 + SIGN(0.5_wp, zQnet) ! Qt > 0 => cooling of the layer => ztf = 1 
    104             !                                 Qt < 0 => warming of the layer => ztf = 0 
    105  
    106             !! Term alpha*Qb (Qb is the virtual surface cooling inc. buoyancy effect of salinity due to evap): 
    107             zz1 = zalpha_w*zQnet - 0.026*MIN(pQlat(ji,jj),0._wp)*rCp0_w/rLevap  ! alpha*(Eq.8) == alpha*Qb "-" because Qlat < 0 
    108             !! LB: this terms only makes sense if > 0 i.e. in the cooling case 
    109             !! so similar to what's done in ECMWF: 
    110             zz1 = MAX(0._wp , zz1)    ! 1. instead of 0.1 though ZQ = MAX(1.0,-pQlw(ji,jj) - pQsen(ji,jj) - pQlat(ji,jj)) 
    111  
    112             zus = MAX(pustar(ji,jj), 1.E-4_wp) ! Laurent: too low wind (u*) might cause problem in stable cases: 
    113             zz2 = zus*zus * roadrw 
    114             zz2 = zz2*zz2 
    115             zlamb =  6._wp*( 1._wp + (rcst_cs*zz1/zz2)**0.75 )**(-1./3.) ! Lambda (Eq.14) (Saunders) 
    116  
    117             ! Updating molecular sublayer thickness (delta): 
    118             zz2    = rnu0_w/(SQRT(roadrw)*zus) 
    119             zdelta =      ztf    *          zlamb*zz2   &  ! Eq.12 (when alpha*Qb>0 / cooling of layer) 
    120                &    + (1._wp - ztf) * MIN(0.007_wp , 6._wp*zz2 )    ! Eq.12 (when alpha*Qb<0 / warming of layer) 
    121             !LB: changed 0.01 to 0.007 
    122  
    123             !! Once again with the new zdelta: 
    124             zfr   = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption / Eq.16 (Fairall al. 1996b) 
    125             zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface 
    126  
    127             !! Update! 
    128             pdT(ji,jj) =  MIN( - zQnet*zdelta/rk0_w , 0._wp )   ! temperature increment 
    129             delta_vl(ji,jj) = zdelta 
     94            zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$ 
     95            !                                       ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$ 
     96 
     97            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 
     98 
     99            DO jc = 1, 4 ! because implicit in terms of zdelta... 
     100               zfr    = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) /  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 ! 
     101               zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface 
     102               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 
     103            END DO 
     104 
     105            dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp )   ! temperature increment 
    130106 
    131107         END DO 
     
    135111 
    136112 
    137  
    138  
    139    SUBROUTINE WL_COARE( pQsw, pQnsol, pTau, pSST, iwait,  pdT, Hwl, mask_wl ) 
     113   SUBROUTINE WL_COARE( pQsw, pQnsol, pTau, pSST, iwait ) 
    140114      !!--------------------------------------------------------------------- 
    141115      !! 
     
    149123      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K] 
    150124      !!     *iwait*      if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop... 
    151       !! 
    152       !!  **   OUTPUT: 
    153       !!     *pdT*   dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST 
    154       !!--------------------------------------------------------------------- 
    155       !! 
    156       !!   ** OPTIONAL OUTPUT: 
    157       !!     *Hwl*        depth of warm layer [m] 
    158       !!     *mask_wl*    mask for possible existence of a warm-layer (1) or not (0) 
    159       !! 
    160       !!------------------------------------------------------------------ 
     125      !!--------------------------------------------------------------------- 
    161126      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 ! 
    162127      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 ! 
     
    164129      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K] 
    165130      INTEGER ,                     INTENT(in)  :: iwait    ! if /= 0 then wait before updating accumulated fluxes 
    166       REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdT      ! dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST 
    167       !! 
    168       REAL(wp),   DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: Hwl     ! depth of warm layer [m] 
    169       INTEGER(1), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: mask_wl ! mask for possible existence of a warm-layer (1) or not (0) 
    170       ! 
    171       ! 
     131      !! 
    172132      INTEGER :: ji,jj 
    173133      ! 
    174       REAL(wp) :: zdT_wl, zQabs, zfr, zdz 
     134      REAL(wp) :: zdTwl, zHwl, zQabs, zfr 
    175135      REAL(wp) :: zqac, ztac 
    176       REAL(wp) :: zalpha_w, zcd1, zcd2, flg 
     136      REAL(wp) :: zalpha, zcd1, zcd2, flg 
    177137      !!--------------------------------------------------------------------- 
    178138 
     
    180140      INTEGER  :: jl 
    181141 
     142      LOGICAL :: l_exit, l_destroy_wl 
     143 
    182144      !! INITIALIZATION: 
    183       pdT(:,:) = 0._wp    ! dT initially set to 0._wp 
    184       zdT_wl = 0._wp       ! total warming (amplitude) in warm layer 
    185       zQabs  = 0._wp       ! total heat absorped in warm layer 
    186       zfr    = zfr0        ! initial value of solar flux absorption 
    187       ztac   = 0._wp 
    188       zqac   = 0._wp 
    189       IF ( PRESENT(mask_wl) ) mask_wl(:,:) = 0 
     145      zQabs  = 0._wp       ! total heat flux absorped in warm layer 
     146      zfr    = zfr0        ! initial value of solar flux absorption !LOLO: save it and use previous value !!! 
    190147 
    191148      IF( .NOT. ln_dm2dc ) CALL sbc_dcy_param() ! we need to call sbc_dcy_param (sbcdcy.F90) because rdawn_dcy and rdusk_dcy are unkonwn otherwize! 
     
    201158         DO ji = 1, jpi 
    202159 
    203             zdz   = MAX( MIN(H_wl(ji,jj),H_wl_max) , 0.1_wp) ! depth of warm layer 
     160            l_exit       = .FALSE. 
     161            l_destroy_wl = .FALSE. 
     162 
     163            zdTwl =  dT_wl(ji,jj)                          ! value of previous time step as first guess 
     164            zHwl  = MAX( MIN(Hz_wl(ji,jj),Hwl_max),0.1_wp) !   "                  "           " 
     165 
     166            zqac = Qnt_ac(ji,jj) ! previous time step Qnt_ac 
     167            ztac = Tau_ac(ji,jj) 
    204168 
    205169            !*****  variables for warm layer  *** 
    206             zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 
    207  
    208             zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha_w*grav*rho0_w))        !mess-o-constants 1 
    209             zcd2 = SQRT(2._wp*zalpha_w*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2 
    210  
    211             !******************************************************** 
    212             !****  Compute apply warm layer  correction ************* 
    213             !******************************************************** 
     170            zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 
     171 
     172            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha*grav*rho0_w))        !mess-o-constants 1 
     173            zcd2 = SQRT(2._wp*zalpha*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2 
     174 
    214175             
    215176            znoon = MOD( 0.5_wp*(rdawn_dcy(ji,jj)+rdusk_dcy(ji,jj)), 1._wp )   ! 0<rnoon<1. => rnoon*24 = UTC time of local noon 
    216177            zmidn = MOD(              znoon-0.5_wp                 , 1._wp ) 
    217  
    218             IF (lwp .AND. (ji==2) .AND. (jj==2) ) WRITE(numout,*) 'LOLO: rdawn, rdusk, znoon, zmidn =', rdawn_dcy(ji,jj), rdusk_dcy(ji,jj), znoon, zmidn 
    219              
    220             ! When midnight is past and dawn is not there yet, do what they call the "midnight reset": 
    221             IF ( (ztime >= zmidn).AND.(ztime <= rdawn_dcy(ji,jj)) ) THEN 
    222                IF (lwp .AND. (ji==2) .AND. (jj==2) ) WRITE(numout,*) ' LOLO [WL_COARE] MIDNIGHT RESET !!!!, ztime =', ztime 
    223                zdz           = H_wl_max 
    224                Tau_ac(ji,jj) = 0._wp 
    225                Qnt_ac(ji,jj) = 0._wp 
    226             END IF 
    227  
    228  
    229             IF ( ztime > rdawn_dcy(ji,jj) ) THEN ! Dawn, a new day starts at current location 
    230                IF (lwp .AND. (ji==2) .AND. (jj==2) ) WRITE(numout,*) ' LOLO [WL_COARE] WE DO WL !!!!, ztime =', ztime 
    231                 
    232                !************************************ 
    233                !****   set warm layer constants  *** 
    234                !************************************ 
    235  
    236                zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)       ! tot heat absorbed in warm layer 
    237  
    238                !PRINT *, '  [WL_COARE] rdt,  pQsw, pQnsol, zQabs =', rdt, REAL(pQsw(ji,jj),4), REAL(pQnsol(ji,jj),4), REAL(zQabs,4) 
    239  
    240                IF ( zQabs >= Qabs_thr ) THEN         ! Check for threshold 
    241  
    242                   !PRINT *, '  [WL_COARE] Tau_ac, Qnt_ac =', REAL(Tau_ac(ji,jj),4), REAL(Qnt_ac(ji,jj),4) 
    243  
    244                   !Tau_ac(ji,jj) = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! momentum integral 
    245                   ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral 
    246  
    247                   IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) THEN         !check threshold for warm layer existence 
    248                      !****************************************** 
    249                      ! Compute the absorption profile 
    250                      !****************************************** 
    251                      DO jl = 1, 5                           !loop 5 times for zfr 
    252                         zfr = 1. - ( 0.28*0.014*(1. - EXP(-zdz/0.014)) + 0.27*0.357*(1. - EXP(-zdz/0.357)) & 
    253                            &        + 0.45*12.82*(1-EXP(-zdz/12.82)) ) / zdz 
    254                         zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed 
    255                         IF ( zqac > 1._wp ) zdz = MAX( MIN( H_wl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth 
    256                      END DO 
    257  
    258                   ELSE 
    259                      !*********************** 
    260                      ! Warm layer wiped out 
    261                      !*********************** 
    262                      zfr  = 0.75 
    263                      zdz  = H_wl_max 
    264                      zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed 
    265  
    266                   END IF !IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) 
    267  
    268                   IF ( zqac > 1._wp ) zdT_wl = zcd2*zqac**1.5/ztac * MAX(zqac/ABS(zqac),0._wp)  !! => IF(zqac>0._wp): zdT_wl=zcd2*zqac**1.5/ztac ; ELSE: zdT_wl=0. / ! normally: zqac > 0 ! 
    269                    
     178            zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight  
     179 
     180            IF ( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN 
     181               ! Dawn reset to 0! 
     182               l_exit       = .TRUE. 
     183               l_destroy_wl = .TRUE. 
     184            END IF 
     185 
     186            IF ( .NOT. l_exit ) THEN 
     187               !! Initial test on initial guess of absorbed heat flux in warm-layer: 
     188               zfr = 1._wp - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & 
     189                  &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 
     190               zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer !LOLO: depends of zfr, which is wild guess... Wrong!!! 
     191               PRINT *, '#LBD:  Initial Qsw & Qnsol:', NINT(pQsw(ji,jj)), NINT(pQnsol(ji,jj)) 
     192               PRINT *, '#LBD:       =>Qabs:', zQabs,' zfr=', zfr 
     193 
     194               IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN 
     195                  ! We have not started to build a WL yet (dT==0) and there's no way it can occur now 
     196                  ! since zQabs <= 0._wp 
     197                  ! => no need to go further 
     198                  PRINT *, '#LBD: we have not started to to build a WL yet (dT==0)' 
     199                  PRINT *, '#LBD: and theres no way it can occur now since zQabs=', zQabs 
     200                  PRINT *, '#LBD: => leaving without changing anything...' 
     201                  l_exit = .TRUE. 
     202               END IF 
     203 
     204            END IF 
     205 
     206            ! Okay test on updated absorbed flux: 
     207            !LOLO: remove??? has a strong influence !!! 
     208            IF ( (.NOT.(l_exit)) .AND. (Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN 
     209               PRINT *, '#LBD: Oh boy! Next Qnt_ac looking weak! =>', Qnt_ac(ji,jj) + zQabs*rdt 
     210               PRINT *, '#LBD:  => time to destroy the warm-layer!' 
     211               l_exit       = .TRUE. 
     212               l_destroy_wl = .TRUE. 
     213            END IF 
     214 
     215 
     216            IF ( .NOT. l_exit) THEN 
     217 
     218               ! Two possibilities at this point: 
     219               ! 1/ A warm layer already exists (dT>0) but it is cooling down because Qabs<0 
     220               ! 2/ Regardless of WL formed (dT==0 or dT>0), we are in the process to initiate one or warm further it ! 
     221 
     222               PRINT *, '#LBD:======================================================' 
     223               PRINT *, '#LBD: WL action makes sense now! => zQabs,dT_wl=', REAL(zQabs,4), REAL(zdTwl,4) 
     224               PRINT *, '#LBD:======================================================' 
     225               PRINT *, '#LBD: current values for Qac and Tac=', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4) 
     226 
     227               ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral 
     228               PRINT *, '#LBD: updated value for Tac=',  REAL(ztac,4) 
     229 
     230               !! We update the value of absorbtion and zQabs: 
     231               !! some part is useless if Qsw=0 !!! 
     232               DO jl = 1, 5 
     233                  zfr = 1. - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & 
     234                     &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 
     235                  zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) 
     236                  zqac  = Qnt_ac(ji,jj) + zQabs*rdt ! updated heat absorbed 
     237                  IF ( zqac <= 0._wp ) EXIT 
     238                  zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth 
     239               END DO 
     240               PRINT *, '#LBD: updated absorption and WL depth=',  REAL(zfr,4), REAL(zHwl,4) 
     241               PRINT *, '#LBD: updated value for Qabs=',  REAL(zQabs,4), 'W/m2' 
     242               PRINT *, '#LBD: updated value for Qac =',  REAL(zqac,4), 'J' 
     243 
     244               IF ( zqac <= 0._wp ) THEN 
     245                  l_destroy_wl = .TRUE. 
     246                  l_exit       = .TRUE. 
     247               ELSE 
     248                  zdTwl = zcd2*zqac**1.5/ztac * MAX(zqac/ABS(zqac),0._wp)  !! => IF(zqac>0._wp): zdTwl=zcd2*zqac**1.5/ztac ; ELSE: zdTwl=0. / ! normally: zqac > 0 ! 
     249                  PRINT *, '#LBD: updated preliminary value for dT_wl=',  REAL(zdTwl,4) 
    270250                  ! Warm layer correction 
    271                   flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = zdT_wl) | 0 when gdept_1d(1)<zdz (pdT(ji,jj) = zdT_wl*gdept_1d(1)/zdz) 
    272                   pdT(ji,jj) = zdT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz ) 
    273                    
    274                END IF ! IF ( zQabs >= Qabs_thr ) 
    275  
    276             END IF ! IF ( isd_sol >= 21600 ) THEN  ! (21600 == 6am) 
     251                  flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl )               ! => 1 when gdept_1d(1)>zHwl (zdTwl = zdTwl) | 0 when gdept_1d(1)<zHwl (zdTwl = zdTwl*gdept_1d(1)/zHwl) 
     252                  zdTwl = zdTwl * ( flg + (1._wp-flg)*gdept_1d(1)/zHwl ) 
     253               END IF 
     254 
     255            END IF !IF ( .NOT. l_exit) 
     256 
     257            IF ( l_destroy_wl ) THEN 
     258               zdTwl = 0._wp 
     259               zfr   = 0.75_wp 
     260               zHwl  = Hwl_max 
     261               zqac  = 0._wp 
     262               ztac  = 0._wp 
     263            END IF 
     264 
     265            PRINT *, '#LBD: exit values for Qac & Tac:', REAL(zqac,4), REAL(ztac,4) 
    277266 
    278267            IF ( iwait == 0 ) THEN 
    279                IF ( (zQabs >= Qabs_thr).AND.(ztime > rdawn_dcy(ji,jj)) ) THEN 
    280                   !PRINT *, '  [WL_COARE] WE UPDATE ACCUMULATED FLUXES !!!' 
    281                   Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral 
    282                   Tau_ac(ji,jj) = ztac ! 
    283                   IF ( PRESENT(mask_wl) ) mask_wl(ji,jj) = 1 
    284                END IF 
    285             END IF 
    286  
    287             H_wl(ji,jj) = zdz 
    288  
    289             IF ( PRESENT(Hwl) ) Hwl(ji,jj) = H_wl(ji,jj) 
     268               !! Iteration loop within bulk algo is over, time to update what needs to be updated: 
     269               dT_wl(ji,jj)  = zdTwl 
     270               Hz_wl(ji,jj)  = zHwl 
     271               PRINT *, '#LBD: FINAL EXIT values for dT_wl & Hz_wl:', REAL(dT_wl(ji,jj),4), REAL(Hz_wl(ji,jj),4) 
     272               Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral 
     273               Tau_ac(ji,jj) = ztac 
     274               PRINT *, '#LBD: FINAL EXIT values for Qac & Tac:', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4) 
     275               PRINT *, '#LBD' 
     276            END IF 
    290277 
    291278         END DO 
    292279      END DO 
    293        
     280 
    294281   END SUBROUTINE WL_COARE 
    295282 
     283 
     284 
     285 
     286   FUNCTION delta_skin_layer( palpha, pQabs, pQlat, pustar_a ) 
     287      !!--------------------------------------------------------------------- 
     288      !! Computes the thickness (m) of the viscous skin layer. 
     289      !! Based on Fairall et al., 1996 
     290      !! 
     291      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 
     292      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 
     293      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 
     294      !! doi:10.1029/95JC03190. 
     295      !! 
     296      !! L. Brodeau, october 2019 
     297      !!--------------------------------------------------------------------- 
     298      REAL(wp)                :: delta_skin_layer 
     299      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!) 
     300      REAL(wp), INTENT(in)    :: pQabs    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 
     301      REAL(wp), INTENT(in)    :: pQlat    ! latent heat flux [W/m^2] 
     302      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s] 
     303      !!--------------------------------------------------------------------- 
     304      REAL(wp) :: zusw, zusw2, zlamb, zQb 
     305      !!--------------------------------------------------------------------- 
     306 
     307      zQb = pQabs + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! LOLO: Double check sign + division by palpha !!! units are okay! 
     308 
     309      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water 
     310      zusw2 = zusw*zusw 
     311 
     312      zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 
     313 
     314      delta_skin_layer = zlamb*rnu0_w/zusw 
     315 
     316   END FUNCTION delta_skin_layer 
    296317 
    297318   !!====================================================================== 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90

    r11637 r11772  
    3838 
    3939   !! Cool-skin related parameters: 
    40    ! ... 
     40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 
     41      &                        dT_cs         !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta) 
     42   REAL(wp), PARAMETER :: zcon0 = -16._wp * grav * rho0_w * rCp0_w * rnu0_w*rnu0_w*rnu0_w / ( rk0_w*rk0_w ) ! "-" because ocean convention: Qabs > 0 => gain of heat for ocean! 
     43   !!                             => see eq.(14) in Fairall et al. 1996   (eq.(6) of Zeng aand Beljaars is WRONG! (typo?) 
    4144 
    4245   !! Warm-layer related parameters: 
    43    REAL(wp), PARAMETER :: rd0  = 3.        !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005) 
    44    REAL(wp), PARAMETER :: zRhoCp_w = rho0_w*rCp0_w 
    45    REAL(wp), PARAMETER :: rNu0 = 1.0       !:  be closer to COARE3p6 ???!LOLO 
    46    !REAL(wp), PARAMETER :: rNu0 = 0.5       !: Nu (exponent of temperature profile) Eq.11 
    47    !                                       !: (Zeng & Beljaars 2005) !: set to 0.5 instead of 
    48    !                                       !: 0.3 to respect a warming of +3 K in calm 
    49    !                                       !: condition for the insolation peak of +1000W/m^2 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 
     47      &                        dT_wl,  &   !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 
     48      &                        Hz_wl       !: depth of warm-layer [m] 
     49   ! 
     50   REAL(wp), PARAMETER, PUBLIC :: rd0  = 3.    !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005) 
     51   REAL(wp), PARAMETER         :: zRhoCp_w = rho0_w*rCp0_w 
     52   ! 
     53   REAL(wp), PARAMETER         :: rNuwl0 = 0.5  !: Nu (exponent of temperature profile) Eq.11 
     54   !                                            !: (Zeng & Beljaars 2005) !: set to 0.5 instead of 
     55   !                                            !: 0.3 to respect a warming of +3 K in calm 
     56   !                                            !: condition for the insolation peak of +1000W/m^2 
    5057   !!---------------------------------------------------------------------- 
    5158CONTAINS 
    5259 
    5360 
    54    SUBROUTINE CS_ECMWF( pQsw, pQnsol, pustar, pSST,  pdT ) 
    55       !!--------------------------------------------------------------------- 
    56       !! 
    57       !!  Cool-Skin scheme according to Fairall et al. 1996 
    58       !! Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL) 
    59       !!  " A prognostic scheme of sea surface skin temperature for modeling and data assimilation " 
    60       !!     as parameterized in IFS Cy45r1 / E.C.M.W.F. 
    61       !!     ------------------------------------------------------------------ 
     61   SUBROUTINE CS_ECMWF( pQsw, pQnsol, pustar, pSST ) 
     62      !!--------------------------------------------------------------------- 
     63      !! 
     64      !! Cool-skin parameterization, based on Fairall et al., 1996: 
     65      !! 
     66      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 
     67      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 
     68      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 
     69      !! doi:10.1029/95JC03190. 
     70      !! 
     71      !!------------------------------------------------------------------ 
    6272      !! 
    6373      !!  **   INPUT: 
     
    6676      !!     *pustar*     friction velocity u*                           [m/s] 
    6777      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K] 
    68       !! 
    69       !!  **  INPUT/OUTPUT: 
    70       !!     *pdT*  : as input =>  previous estimate of dT cool-skin 
    71       !!              as output =>  new estimate of dT cool-skin 
    72       !! 
    7378      !!------------------------------------------------------------------ 
    74       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 
    75       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 
    76       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*) 
    77       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pSST ! bulk SST [K] 
    78       !! 
    79       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT    ! dT due to cool-skin effect 
    80       !!--------------------------------------------------------------------- 
    81       INTEGER  ::   ji, jj     ! dummy loop indices 
    82       REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, & 
    83          & zusw, zusw2 
    84       !!--------------------------------------------------------------------- 
    85  
     79      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 
     80      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 
     81      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*) 
     82      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K] 
     83      !!--------------------------------------------------------------------- 
     84      INTEGER  :: ji, jj, jc 
     85      REAL(wp) :: zQabs, zdelta, zfr 
     86      !!--------------------------------------------------------------------- 
    8687      DO jj = 1, jpj 
    8788         DO ji = 1, jpi 
    8889 
    89             zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 
    90  
    91             zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere 
    92  
    93             zusw  = MAX(pustar(ji,jj), 1.E-4_wp)*SQRT(roadrw)    ! u* in the water 
    94             zusw2 = zusw*zusw 
    95  
    96             zlamb = 6._wp*( 1._wp + (zQnsol*zalpha_w*rcst_cs/(zusw2*zusw2 ))**0.75 )**(-1./3.)   ! w.r.t COARE 3p6 => seems to ommit absorbed zfr*Qsw (Qnet i.o. Qnsol) and effect of evap... 
    97             !                                                                                  ! so zlamb not implicit in terms of zdelta (zfr(delta)), so no need to have last guess of delta as in COARE 3.6... 
    98             zdelta = zlamb*rnu0_w/zusw 
    99  
    100             zfr   = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption / Eq. 8.131 / IFS cy40r1, doc, Part IV, 
    101             zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface 
    102  
    103             !! Update! 
    104             pdT(ji,jj) =  MIN( - zQnet*zdelta/rk0_w , 0._wp )   ! temperature increment 
     90            zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$ 
     91            !                                       ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$ 
     92            !zQabs  = pQnsol(ji,jj) 
     93 
     94            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 
     95 
     96            DO jc = 1, 4 ! because implicit in terms of zdelta... 
     97               zfr    = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005 
     98               !              =>  (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996) 
     99               zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface 
     100               !zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) 
     101               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 
     102            END DO 
     103 
     104            dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp )   ! temperature increment 
    105105 
    106106         END DO 
     
    111111 
    112112 
    113    SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST, pdT ) 
     113   SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST,  pustk ) 
    114114      !!--------------------------------------------------------------------- 
    115115      !! 
    116116      !!  Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL) 
    117117      !!  " A prognostic scheme of sea surface skin temperature for modeling and data assimilation " 
     118      !! 
     119      !!  STIL NO PROGNOSTIC EQUATION FOR THE DEPTH OF THE WARM-LAYER! 
    118120      !! 
    119121      !!    As included in IFS Cy45r1   /  E.C.M.W.F. 
     
    125127      !!     *pustar*     friction velocity u*                           [m/s] 
    126128      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K] 
    127       !! 
    128       !!   **  INPUT/OUTPUT: 
    129       !!     *pdT*  : as input =>  previous estimate of dT warm-layer 
    130       !!             as output =>  new estimate of dT warm-layer 
    131       !! 
    132129      !!------------------------------------------------------------------ 
    133130      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 ! 
     
    135132      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pustar   ! friction velocity [m/s] 
    136133      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K] 
    137       ! 
    138       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT    ! dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST 
    139       ! 
    140       INTEGER :: ji,jj 
     134      !! 
     135      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pustk ! surface Stokes velocity [m/s] 
     136      ! 
     137      INTEGER :: ji, jj, jc 
    141138      ! 
    142139      REAL(wp) :: & 
    143          & zdz,    & !: thickness of the warm-layer [m] 
    144          & zalpha_w, & !: thermal expansion coefficient of sea-water 
    145          & ZSRD,   & 
    146          & dT_wl,   & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL 
    147          & zfr,zdL,zdL2, ztmp, & 
    148          & ZPHI, & 
    149          & zus_a, zusw, zusw2, & 
    150          & flg, zQabs, ZL1, ZL2 
    151       !!--------------------------------------------------------------------- 
     140         & zHwl,    &  !: thickness of the warm-layer [m] 
     141         & ztcorr,  &  !: correction of dT w.r.t measurement depth of bulk SST (first T-point) 
     142         & zalpha, & !: thermal expansion coefficient of sea-water [1/K] 
     143         & zdTwl_b, zdTwl_n, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL 
     144         & zfr, zeta, ztmp, & 
     145         & zusw, zusw2, & 
     146         & zLa, zfLa, & 
     147         & flg, zwf, zQabs, & 
     148         & ZA, ZB, zL1, zL2, & 
     149         &  zcst0, zcst1, zcst2, zcst3 
     150      ! 
     151      LOGICAL :: l_pustk_known 
     152      !!--------------------------------------------------------------------- 
     153 
     154      l_pustk_known = .FALSE. 
     155      IF ( PRESENT(pustk) ) l_pustk_known = .TRUE. 
    152156 
    153157      DO jj = 1, jpj 
    154158         DO ji = 1, jpi 
    155159 
    156             zdz = rd0 ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !) 
    157  
    158             ! dT_wl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zdz) 
     160            zHwl = Hz_wl(ji,jj) ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !) 
     161            ! it is = rd0 (3m) in default Zeng & Beljaars case... 
     162 
     163            !! Previous value of dT / warm-layer, adapted to depth: 
     164            flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl )               ! => 1 when gdept_1d(1)>zHwl (dT_wl(ji,jj) = zdTwl) | 0 when z_s$ 
     165            ztcorr = flg + (1._wp - flg)*gdept_1d(1)/zHwl 
     166            zdTwl_b = MAX ( dT_wl(ji,jj) / ztcorr , 0._wp ) 
     167            ! zdTwl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zHwl) 
    159168            ! pdT         "                          "                                    and depth of bulk SST (here gdept_1d(1))! 
    160             !! => but of course in general the bulk SST is taken shallower than zdz !!! So correction less pronounced! 
    161             !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof dT_wl ! 
    162             flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = dT_wl) | 0 when z_s$ 
    163             dT_wl = pdT(ji,jj) / ( flg + (1._wp-flg)*gdept_1d(1)/zdz ) 
    164  
    165             zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 
     169            !! => but of course in general the bulk SST is taken shallower than zHwl !!! So correction less pronounced! 
     170            !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof zdTwl ! 
     171 
     172            zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 
    166173 
    167174 
    168175            ! *** zfr = Fraction of solar radiation absorbed in warm layer (-) 
    169             zfr = 1._wp - 0.28_wp*EXP(-71.5_wp*zdz) - 0.27_wp*EXP(-2.8_wp*zdz) - 0.45_wp*EXP(-0.07_wp*zdz)  !: Eq. 8.157 
     176            zfr = 1._wp - 0.28_wp*EXP(-71.5_wp*zHwl) - 0.27_wp*EXP(-2.8_wp*zHwl) - 0.45_wp*EXP(-0.07_wp*zHwl)  !: Eq. 8.157 
    170177 
    171178            zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)       ! tot heat absorbed in warm layer 
    172179 
    173             zusw  = MAX(pustar(ji,jj), 1.E-4_wp)*SQRT(roadrw)    ! u* in the water 
     180            zusw  = MAX( pustar(ji,jj), 1.E-4_wp ) * sq_radrw    ! u* in the water 
    174181            zusw2 = zusw*zusw 
    175182 
    176  
    177             !! *** 1st rhs term in eq. 8.156 (IFS doc Cy45r1): 
    178             ZL1 = zQabs / ( zdz * zRhoCp_w * rNu0 ) * (rNu0 + 1._wp) 
    179  
    180  
    181             !! Buoyancy flux and stability parameter (zdl = -z/L) in water 
    182             ZSRD = zQabs/zRhoCp_w 
    183             ! 
    184             flg = 0.5_wp + SIGN(0.5_wp, ZSRD)  ! ZSRD > 0. => 1.  / ZSRD < 0. => 0. 
    185             ztmp = MAX(dT_wl,0._wp) 
    186             zdl = (flg+1._wp) * ( zusw2 * SQRT(ztmp/(5._wp*zdz*grav*zalpha_w/rNu0)) ) & ! (dT_wl > 0.0 .AND. ZSRD < 0.0) 
    187                &  +    flg    *  ZSRD                                                                  !   otherwize 
    188             ! 
    189             zus_a = MAX( pustar(ji,jj), 1.E-4_wp ) 
    190             zdL = zdz*vkarmn*grav/(roadrw)**1.5_wp*zalpha_w*zdL/(zus_a*zus_a*zus_a) 
    191  
    192             !! Stability function Phi_t(-z/L) (zdL is -z/L) : 
    193             flg = 0.5_wp + SIGN(0.5_wp, zdL)  ! zdl > 0. => 1.  / zdl < 0. => 0. 
    194             zdL2 = zdL*zdL 
    195             ZPHI =    flg      * ( 1._wp + (5._wp*zdL + 4._wp*zdL2)/(1._wp + 3._wp*zdL + 0.25_wp*zdL2) ) &  ! (zdL > 0) Takaya et al. 
    196                & + (flg+1._wp) * ( 1._wp/SQRT(1._wp - 16._wp*(-ABS(zdL))) )        ! (zdl < 0) Eq. 8.136 
    197             !! FOR zdL > 0.0, old relations: 
    198             !         ZPHI = 1.+5._wp*zdL                                ! Eq. 8.136 (Large et al. 1994) 
    199             !         ZPHI = 1.+5.0*(zdL+zdL**2)/(1.0+3.0*zdL+zdL**2) ! SHEBA, Grachev et al. 2007 
    200  
    201             !! *** 2nd rhs term in eq. 8.156 (IFS doc Cy45r1): 
    202             ZL2 = - (rNu0 + 1._wp) * vkarmn * zusw / ( zdz * ZPHI ) 
    203  
    204             ! Forward time / explicit solving of eq. 8.156 (IFS doc Cy45r1): (f_n+1 == pdT(ji,jj) ; f_n == dT_wl) 
    205             dT_wl = MAX ( dT_wl + rdt*ZL1 + rdt*ZL2*dT_wl , 0._wp ) 
    206  
    207             ! dT_wl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zdz) 
    208             !! => but of course in general the bulk SST is taken shallower than zdz !!! So correction less pronounced! 
    209  
    210             flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = dT_wl) | 0 when z_s$ 
    211             pdT(ji,jj) = dT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz ) 
    212  
     183            ! Langmuir: 
     184            IF ( l_pustk_known ) THEN 
     185               zLa = SQRT(zusw/MAX(pustk(ji,jj),1.E-6)) 
     186            ELSE 
     187               zla = 0.3_wp 
     188            END IF 
     189            zfLa = MAX( zla**(-2._wp/3._wp) , 1._wp )   ! Eq.(6) 
     190 
     191            zwf = 0.5_wp + SIGN(0.5_wp, zQabs)  ! zQabs > 0. => 1.  / zQabs < 0. => 0. 
     192 
     193            zcst1 = vkarmn*grav*zalpha 
     194 
     195            ! 1/L when zQabs > 0 : 
     196            zL2 = zcst1*zQabs / (zRhoCp_w*zusw2*zusw) 
     197                
     198            zcst2 = zcst1 / ( 5._wp*zHwl*zusw2 )  !OR: zcst2 = zcst1*rNuwl0 / ( 5._wp*zHwl*zusw2 ) ??? 
     199 
     200            zcst0 = rdt * (rNuwl0 + 1._wp) / zHwl 
     201             
     202            ZA = zcst0 * zQabs / ( rNuwl0 * zRhoCp_w ) 
     203 
     204            zcst3 = -zcst0 * vkarmn * zusw * zfLa 
     205 
     206            !! Sorry about all these constants ( constant w.r.t zdTwl), it's for 
     207            !! the sake of optimizations... So all these operations are not done 
     208            !! over and over within the iteration loop... 
     209             
     210            !! T R U L L Y   I M P L I C I T => needs itteration 
     211            !! => have to itterate just because the 1/(Monin-Obukhov length), zL1, uses zdTwl when zQabs < 0.. 
     212            !!    (without this term otherwize the implicit analytical solution is straightforward...) 
     213            zdTwl_n = zdTwl_b 
     214            DO jc = 1, 10 
     215                
     216               zdTwl_n = 0.5_wp * ( zdTwl_n + zdTwl_b ) ! semi implicit, for faster convergence 
     217                
     218               ! 1/L when zdTwl > 0 .AND. zQabs < 0 : 
     219               zL1 =         SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalpha/( 5._wp*zHwl ) ) / zusw 
     220               !zL1 = vkarmn*SQRT( zdTwl_n       *grav*zalpha        / ( 5._wp*zHwl ) ) / zusw   ! => vkarmn outside, not inside zcst1 (just for this particular line) ??? 
     221                
     222               ! Stability parameter (z/L): 
     223               zeta =  (1._wp - zwf) * zHwl*zL1   +   zwf * zHwl*zL2 
     224 
     225               ZB = zcst3 / PHI(zeta) 
     226 
     227               zdTwl_n = MAX ( zdTwl_b + ZA + ZB*zdTwl_n , 0._wp )  ! Eq.(6) 
     228 
     229            END DO 
     230             
     231            !! Update: 
     232            dT_wl(ji,jj) = zdTwl_n * ztcorr 
     233             
    213234         END DO 
    214235      END DO 
     
    216237   END SUBROUTINE WL_ECMWF 
    217238 
     239 
     240 
     241   FUNCTION delta_skin_layer( palpha, pQabs, pustar_a ) 
     242      !!--------------------------------------------------------------------- 
     243      !! Computes the thickness (m) of the viscous skin layer. 
     244      !! Based on Fairall et al., 1996 
     245      !! 
     246      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 
     247      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 
     248      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 
     249      !! doi:10.1029/95JC03190. 
     250      !! 
     251      !! L. Brodeau, october 2019 
     252      !!--------------------------------------------------------------------- 
     253      REAL(wp)                :: delta_skin_layer 
     254      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!) 
     255      REAL(wp), INTENT(in)    :: pQabs    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 
     256      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s] 
     257      !!--------------------------------------------------------------------- 
     258      REAL(wp) :: zusw, zusw2, zlamb, zQb 
     259      !!--------------------------------------------------------------------- 
     260 
     261      zQb = pQabs 
     262 
     263      !zQ = MIN( -0.1_wp , pQabs ) 
     264 
     265      !ztf = 0.5_wp + SIGN(0.5_wp, zQ)  ! Qabs < 0 => cooling of the layer => ztf = 0 (normal case) 
     266      !                                   ! Qabs > 0 => warming of the layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux) 
     267      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water 
     268      zusw2 = zusw*zusw 
     269 
     270      zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 
     271      !zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*zQ, 0._wp)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 
     272 
     273      delta_skin_layer = zlamb*rnu0_w/zusw 
     274 
     275      !delta_skin_layer =  (1._wp - ztf) * zlamb*rnu0_w/zusw    &         ! see eq.(12) in Fairall et al., 1996 
     276      !   &               +     ztf  * MIN(6._wp*rnu0_w/zusw , 0.007_wp) 
     277   END FUNCTION delta_skin_layer 
     278 
     279 
     280 
     281   FUNCTION PHI( pzeta) 
     282      !!--------------------------------------------------------------------- 
     283      !! 
     284      !! Takaya et al., 2010 
     285      !!  Eq.(5) 
     286      !! L. Brodeau, october 2019 
     287      !!--------------------------------------------------------------------- 
     288      REAL(wp)                :: PHI 
     289      REAL(wp), INTENT(in)    :: pzeta    ! stability parameter 
     290      !!--------------------------------------------------------------------- 
     291      REAL(wp) :: ztf, zzt2 
     292      !!--------------------------------------------------------------------- 
     293      ! 
     294      zzt2 = pzeta*pzeta 
     295      ! 
     296      ztf = 0.5_wp + SIGN(0.5_wp, pzeta)  ! zeta > 0 => ztf = 1 
     297      !                                   ! zeta < 0 => ztf = 0 
     298      PHI =      ztf     * ( 1. + (5.*pzeta + 4.*zzt2)/(1. + 3.*pzeta + 0.25*zzt2) ) &   ! zeta > 0 
     299         &  + (1. - ztf) * 1./SQRT( 1. - 16.*(-ABS(pzeta)) )                             ! zeta < 0 
     300      ! 
     301   END FUNCTION PHI 
     302 
    218303   !!====================================================================== 
    219304END MODULE sbcblk_skin_ecmwf 
Note: See TracChangeset for help on using the changeset viewer.