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 11672 – NEMO

Changeset 11672


Ignore:
Timestamp:
2019-10-10T12:32:53+02:00 (3 years ago)
Author:
laurent
Message:

LB: "sbcblk_skin_coare.F90" now relies on "sbcdcy.F90" (modified) to know dawn/dusk time

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

Legend:

Unmodified
Added
Removed

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

    r11665 r11672  
    529529 
    530530         CASE( np_COARE_3p0 ) 
    531             IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p0" WITH CSWL options!!!, nsec_day, gdept_1d(1)=', nsec_day, gdept_1d(1) !LBrm 
     531            IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p0" WITH CSWL options!!!, gdept_1d(1)=', gdept_1d(1) !LBrm 
    532532            CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.0 
    533533               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    & 
    534                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), & 
    535                &                isecday_utc=nsec_day, plong=glamt(:,:) ) 
     534               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
    536535 
    537536         CASE( np_COARE_3p6 ) 
    538             IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p6" WITH CSWL options!!!, nsec_day, gdept_1d(1)=', nsec_day, gdept_1d(1) !LBrm 
     537            IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p6" WITH CSWL options!!!, gdept_1d(1)=', gdept_1d(1) !LBrm 
    539538            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6 
    540539               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    & 
    541                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), & 
    542                &                isecday_utc=nsec_day, plong=glamt(:,:) ) 
     540               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
    543541 
    544542         CASE( np_ECMWF     ) 
    545             IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_ecmwf" WITH CSWL options!!!, nsec_day, gdept_1d(1)=', nsec_day, gdept_1d(1) !LBrm 
     543            IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_ecmwf" WITH CSWL options!!!, gdept_1d(1)=', gdept_1d(1) !LBrm 
    546544            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF 
    547545               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 

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

    r11666 r11672  
    8787 
    8888 
    89    SUBROUTINE turb_coare3p0( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl,  & 
    90       &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                               & 
    91       &                      Cdn, Chn, Cen,                      & 
    92       &                      Qsw, rad_lw, slp, pdT_cs,                                    & ! optionals for cool-skin (and warm-layer) 
    93       &                      isecday_utc, plong, pdT_wl, Hwl )                             ! optionals for warm-layer only 
     89   SUBROUTINE turb_coare3p0( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 
     90      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                              & 
     91      &                      Cdn, Chn, Cen,                                              & 
     92      &                      Qsw, rad_lw, slp, pdT_cs,                                   & ! optionals for cool-skin (and warm-layer) 
     93      &                      pdT_wl, Hwl )                                                 ! optionals for warm-layer only 
    9494      !!---------------------------------------------------------------------- 
    9595      !!                      ***  ROUTINE  turb_coare3p0  *** 
    … …  
    132132      !!    *  slp    : sea-level pressure                                    [Pa] 
    133133      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K] 
    134       !!    * isecday_utc: 
    135       !!    *  plong  : longitude array                                       [deg.E] 
    136134      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K] 
    137135      !!    * Hwl     : depth of warm layer                                   [m] 
    … …  
    172170      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs 
    173171      ! 
    174       INTEGER,  INTENT(in   ), OPTIONAL                     ::   isecday_utc ! current UTC time, counted in second since 00h of the current day 
    175       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   plong    !             [deg.E] 
    176172      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K] 
    177173      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Hwl      !             [m] 
    … …  
    194190         &                pdTw,   &  ! SST increment "dT" for warm layer correction          [K] 
    195191         &                zHwl       ! depth of warm-layer [m] 
    196  
    197       ! 
    198       LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE. 
    199       CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0.F90' 
    200       CHARACTER(len=128) :: cf_tmp 
    201192      !!---------------------------------------------------------------------------------- 
    202193 
    203194      IF ( kt == 1 ) CALL COARE3P0_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays 
    204  
    205195 
    206196      l_zt_equal_zu = .FALSE. 
    … …  
    209199 
    210200      !! Initializations for cool skin and warm layer: 
    211       IF ( l_use_cs ) 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 cool-skin param!' 
    214             STOP 
    215          END IF 
    216          ALLOCATE ( pdTc(jpi,jpj) ) 
    217          pdTc(:,:) = -0.25_wp  ! First guess of skin correction 
    218       END IF 
    219  
    220       IF ( l_use_wl ) THEN 
    221          IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) .AND. PRESENT(isecday_utc) .AND. PRESENT(plong))) THEN 
    222             PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw, slp, isecday_utc & plong to use warm-layer param!' 
    223             STOP 
    224          END IF 
    225          ALLOCATE ( pdTw(jpi,jpj) ) 
    226          IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) ) 
    227       END IF 
    228  
    229201      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' ) 
    230204         ALLOCATE ( zsst(jpi,jpj) ) 
    231205         zsst = T_s ! backing up the bulk SST 
    232206         IF( l_use_cs ) T_s = T_s - 0.25   ! First guess of correction 
    233207         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) ) 
    234216      END IF 
    235217 
    … …  
    358340            !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 
    359341            IF (PRESENT(Hwl)) THEN 
    360                CALL WL_COARE( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl ) 
     342               CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl ) 
    361343            ELSE 
    362                CALL WL_COARE( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw ) 
     344               CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt),  pdTw ) 
    363345            END IF 
    364346            !! Updating T_s and q_s !!! 

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

    r11666 r11672  
    8787 
    8888 
    89    SUBROUTINE turb_coare3p6( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl,  & 
    90       &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                               & 
    91       &                      Cdn, Chn, Cen,                      & 
    92       &                      Qsw, rad_lw, slp, pdT_cs,                                    & ! optionals for cool-skin (and warm-layer) 
    93       &                      isecday_utc, plong, pdT_wl, Hwl )                             ! optionals for warm-layer only 
     89   SUBROUTINE turb_coare3p6( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 
     90      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                              & 
     91      &                      Cdn, Chn, Cen,                                              & 
     92      &                      Qsw, rad_lw, slp, pdT_cs,                                   & ! optionals for cool-skin (and warm-layer) 
     93      &                      pdT_wl, Hwl )                                                 ! optionals for warm-layer only 
    9494      !!---------------------------------------------------------------------- 
    9595      !!                      ***  ROUTINE  turb_coare3p6  *** 
    … …  
    123123      !! 
    124124      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg] 
    125       !!              -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True) 
    126       !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False) 
     125      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True) 
     126      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False) 
    127127      !! 
    128128      !! OPTIONAL INPUT: 
    … …  
    132132      !!    *  slp    : sea-level pressure                                    [Pa] 
    133133      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K] 
    134       !!    * isecday_utc: 
    135       !!    *  plong  : longitude array                                       [deg.E] 
    136134      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K] 
    137135      !!    * Hwl     : depth of warm layer                                   [m] 
    … …  
    172170      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs 
    173171      ! 
    174       INTEGER,  INTENT(in   ), OPTIONAL                     ::   isecday_utc ! current UTC time, counted in second since 00h of the current day 
    175       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   plong    !             [deg.E] 
    176172      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K] 
    177173      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Hwl      !             [m] 
    … …  
    194190         &                pdTw,   &  ! SST increment "dT" for warm layer correction          [K] 
    195191         &                zHwl       ! depth of warm-layer [m] 
    196  
    197       ! 
    198       LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE. 
    199       CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p6@sbcblk_algo_coare3p6.F90' 
    200       CHARACTER(len=128) :: cf_tmp 
    201192      !!---------------------------------------------------------------------------------- 
    202193 
    203194      IF ( kt == 1 ) CALL COARE3P6_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays 
    204  
    205195 
    206196      l_zt_equal_zu = .FALSE. 
    … …  
    209199 
    210200      !! Initializations for cool skin and warm layer: 
    211       IF ( l_use_cs ) 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 cool-skin param!' 
    214             STOP 
    215          END IF 
    216          ALLOCATE ( pdTc(jpi,jpj) ) 
    217          pdTc(:,:) = -0.25_wp  ! First guess of skin correction 
    218       END IF 
    219  
    220       IF ( l_use_wl ) THEN 
    221          IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) .AND. PRESENT(isecday_utc) .AND. PRESENT(plong))) THEN 
    222             PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw, slp, isecday_utc & plong to use warm-layer param!' 
    223             STOP 
    224          END IF 
    225          ALLOCATE ( pdTw(jpi,jpj) ) 
    226          IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) ) 
    227       END IF 
    228  
    229201      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' ) 
    230204         ALLOCATE ( zsst(jpi,jpj) ) 
    231205         zsst = T_s ! backing up the bulk SST 
    232206         IF( l_use_cs ) T_s = T_s - 0.25   ! First guess of correction 
    233207         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) ) 
    234216      END IF 
    235217 
    … …  
    358340            !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 
    359341            IF (PRESENT(Hwl)) THEN 
    360                CALL WL_COARE( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl ) 
     342               CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl ) 
    361343            ELSE 
    362                CALL WL_COARE( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw ) 
     344               CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt),  pdTw ) 
    363345            END IF 
    364346            !! Updating T_s and q_s !!! 

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

    r11626 r11672  
    2626   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc... 
    2727 
     28   USE sbcdcy          !LOLO: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE) 
     29    
    2830   USE lib_mpp         ! distribued memory computing library 
    2931   USE in_out_manager  ! I/O manager 
    … …  
    135137 
    136138 
    137    SUBROUTINE WL_COARE( kt,  pQsw, pQnsol, pTau, pSST, plon, isd, iwait,  pdT, & 
    138       &                    Hwl, mask_wl ) 
     139   SUBROUTINE WL_COARE( pQsw, pQnsol, pTau, pSST, iwait,  pdT, Hwl, mask_wl ) 
    139140      !!--------------------------------------------------------------------- 
    140141      !! 
    … …  
    147148      !!     *pTau*       surface wind stress                            [N/m^2] 
    148149      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K] 
    149       !!     *plon*       longitude                                      [deg.E] 
    150       !!     *isd*        current UTC time, counted in second since 00h of the current day 
    151150      !!     *iwait*      if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop... 
    152151      !! 
    … …  
    160159      !! 
    161160      !!------------------------------------------------------------------ 
    162       INTEGER ,                     INTENT(in)  :: kt 
    163161      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 ! 
    164162      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 ! 
    165163      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTau     ! wind stress [N/m^2] 
    166164      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K] 
    167       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: plon     ! longitude [deg.E] 
    168       INTEGER ,                     INTENT(in)  :: isd      ! current UTC time, counted in second since 00h of the current day 
    169165      INTEGER ,                     INTENT(in)  :: iwait    ! if /= 0 then wait before updating accumulated fluxes 
    170166      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 
    … …  
    179175      REAL(wp) :: zqac, ztac 
    180176      REAL(wp) :: zalpha_w, zcd1, zcd2, flg 
    181       CHARACTER(len=128) :: cf_tmp 
    182       !!--------------------------------------------------------------------- 
    183  
    184       REAL(wp) :: rlag_gw_h, &  ! local solar time lag in hours   / Greenwich meridian (lon==0) => ex: ~ -10.47 hours for Hawai 
    185          &        rhr_sol       ! local solar time in hours since midnight 
    186  
    187       INTEGER  :: ilag_gw_s, &  ! local solar time LAG in seconds / Greenwich meridian (lon==0) => ex: ~ INT( -10.47*3600. ) seconds for Hawai 
    188          &        isd_sol,   &  ! local solar time in number of seconds since local solar midnight 
    189          &        jl 
     177      !!--------------------------------------------------------------------- 
     178 
     179      REAL(wp) :: ztime, znoon, zmidn 
     180      INTEGER  :: jl 
    190181 
    191182      !! INITIALIZATION: 
    … …  
    198189      IF ( PRESENT(mask_wl) ) mask_wl(:,:) = 0 
    199190 
     191      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! 
     192 
     193      ztime = REAL(nsec_day,wp)/(24._wp*3600._wp) ! time of current time step since 00:00 for current day (UTC) -> ztime = 0 -> 00:00 / ztime = 0.5 -> 12:00 ... 
     194 
     195      IF (lwp) THEN 
     196         WRITE(numout,*) '' 
     197         WRITE(numout,*) 'LOLO: sbcblk_skin_coare => nsec_day, ztime =', nsec_day, ztime 
     198      END IF 
     199       
    200200      DO jj = 1, jpj 
    201201         DO ji = 1, jpi 
    202202 
    203203            zdz   = MAX( MIN(H_wl(ji,jj),H_wl_max) , 0.1_wp) ! depth of warm layer 
    204  
    205             !! Need to know the local solar time from longitude and isd: 
    206             rlag_gw_h = -1._wp * MODULO( ( 360._wp - MODULO(plon(ji,jj),360._wp) ) / 15._wp , 24._wp ) 
    207             rlag_gw_h = -1._wp * SIGN( MIN(ABS(rlag_gw_h) , ABS(MODULO(rlag_gw_h,24._wp))), rlag_gw_h + 12._wp ) 
    208             ilag_gw_s = INT( rlag_gw_h*3600._wp ) 
    209             isd_sol = MODULO( isd + ilag_gw_s , 24*3600 ) 
    210             rhr_sol = REAL( isd_sol , wp) / 3600._wp 
    211204 
    212205            !*****  variables for warm layer  *** 
    … …  
    220213            !******************************************************** 
    221214             
    222             !IF (isd_sol <= rdt ) THEN    !re-zero at midnight ! LOLO improve: risky if real midnight (00:00:00) is not a time in vtime... 
    223             IF ( (rhr_sol > 23.5_wp).OR.(rhr_sol < 4._wp) ) THEN 
    224                !PRINT *, '  [WL_COARE] MIDNIGHT RESET !!!!, isd_sol =>', isd_sol 
     215            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 
     216            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 
    225223               zdz           = H_wl_max 
    226224               Tau_ac(ji,jj) = 0._wp 
    … …  
    228226            END IF 
    229227 
    230             IF ( rhr_sol > 5._wp ) THEN  ! ( 5am) 
    231  
    232                !PRINT *, '  [WL_COARE] WE DO WL !' 
    233                !PRINT *, '  [WL_COARE] isd_sol, pTau, pSST, pdT =', isd_sol, REAL(pTau(ji,jj),4), REAL(pSST(ji,jj),4), REAL(pdT(ji,jj),4) 
    234  
     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                
    235232               !************************************ 
    236233               !****   set warm layer constants  *** 
    … …  
    280277 
    281278            IF ( iwait == 0 ) THEN 
    282                IF ( (zQabs >= Qabs_thr).AND.(rhr_sol >= 5._wp) ) THEN 
     279               IF ( (zQabs >= Qabs_thr).AND.(ztime > rdawn_dcy(ji,jj)) ) THEN 
    283280                  !PRINT *, '  [WL_COARE] WE UPDATE ACCUMULATED FLUXES !!!' 
    284281                  Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral 

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

    r10425 r11672  
    77   !!   NEMO    2.0  !  2006-02  (S. Masson, G. Madec)  adaptation to NEMO 
    88   !!           3.1  !  2009-07  (J.M. Molines)  adaptation to v3.1 
     9   !!           4.*  !  2019-10  (L. Brodeau)  nothing really new, but the routine 
     10   !!                ! "sbc_dcy_param" has been extracted from old function "sbc_dcy" 
     11   !!                ! => this allows the warm-layer param of COARE3* to know the time 
     12   !!                ! of dawn and dusk even if "ln_dm2dc=.false." (rdawn_dcy & rdusk_dcy 
     13   !!                ! are now public) 
    914   !!---------------------------------------------------------------------- 
    1015 
    … …  
    2227   IMPLICIT NONE 
    2328   PRIVATE 
    24     
     29 
    2530   INTEGER, PUBLIC ::   nday_qsr   !: day when parameters were computed 
    26     
     31 
    2732   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   raa , rbb  , rcc  , rab     ! diurnal cycle parameters 
    28    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rtmd, rdawn, rdusk, rscal   !    -      -       - 
    29    
     33   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rtmd, rscal   !    -      -       - 
     34   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: rdawn_dcy, rdusk_dcy   !    -      -       - 
     35 
    3036   PUBLIC   sbc_dcy        ! routine called by sbc 
     37   PUBLIC   sbc_dcy_param  ! routine used here and called by warm-layer parameterization (sbcblk_skin_coare*) 
    3138 
    3239   !!---------------------------------------------------------------------- 
    3340   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    34    !! $Id$  
     41   !! $Id$ 
    3542   !! Software governed by the CeCILL license (see ./LICENSE) 
    3643   !!---------------------------------------------------------------------- 
    3744CONTAINS 
    3845 
    39       INTEGER FUNCTION sbc_dcy_alloc() 
    40          !!---------------------------------------------------------------------- 
    41          !!                ***  FUNCTION sbc_dcy_alloc  *** 
    42          !!---------------------------------------------------------------------- 
    43          ALLOCATE( raa (jpi,jpj) , rbb  (jpi,jpj) , rcc  (jpi,jpj) , rab  (jpi,jpj) ,     & 
    44             &      rtmd(jpi,jpj) , rdawn(jpi,jpj) , rdusk(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc ) 
    45             ! 
    46          CALL mpp_sum ( 'sbcdcy', sbc_dcy_alloc ) 
    47          IF( sbc_dcy_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc: failed to allocate arrays' ) 
    48       END FUNCTION sbc_dcy_alloc 
     46   INTEGER FUNCTION sbc_dcy_alloc() 
     47      !!---------------------------------------------------------------------- 
     48      !!                ***  FUNCTION sbc_dcy_alloc  *** 
     49      !!---------------------------------------------------------------------- 
     50      ALLOCATE( raa (jpi,jpj) , rbb  (jpi,jpj) , rcc  (jpi,jpj) , rab  (jpi,jpj) ,     & 
     51         &      rtmd(jpi,jpj) , rdawn_dcy(jpi,jpj) , rdusk_dcy(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc ) 
     52      ! 
     53      CALL mpp_sum ( 'sbcdcy', sbc_dcy_alloc ) 
     54      IF( sbc_dcy_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc: failed to allocate arrays' ) 
     55   END FUNCTION sbc_dcy_alloc 
    4956 
    5057 
    … …  
    6067      !! 
    6168      !! reference  : Bernie, DJ, E Guilyardi, G Madec, JM Slingo, and SJ Woolnough, 2007 
    62       !!              Impact of resolving the diurnal cycle in an ocean--atmosphere GCM.  
     69      !!              Impact of resolving the diurnal cycle in an ocean--atmosphere GCM. 
    6370      !!              Part 1: a diurnally forced OGCM. Climate Dynamics 29:6, 575-590. 
    6471      !!---------------------------------------------------------------------- 
    6572      LOGICAL , OPTIONAL          , INTENT(in) ::   l_mask    ! use the routine for night mask computation 
    66       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqsrin    ! input daily QSR flux  
     73      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqsrin    ! input daily QSR flux 
    6774      REAL(wp), DIMENSION(jpi,jpj)             ::   zqsrout   ! output QSR flux with diurnal cycle 
    6875      !! 
    6976      INTEGER  ::   ji, jj                                       ! dummy loop indices 
    7077      INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask 
    71       REAL(wp) ::   ztwopi, zinvtwopi, zconvrad  
    7278      REAL(wp) ::   zlo, zup, zlousd, zupusd 
    73       REAL(wp) ::   zdsws, zdecrad, ztx, zsin, zcos 
    74       REAL(wp) ::   ztmp, ztmp1, ztmp2, ztest 
     79      REAL(wp) ::   ztmp, ztmp1, ztmp2 
    7580      REAL(wp) ::   ztmpm, ztmpm1, ztmpm2 
    76       !---------------------------statement functions------------------------ 
    77       REAL(wp) ::   fintegral, pt1, pt2, paaa, pbbb, pccc        ! dummy statement function arguments 
    78       fintegral( pt1, pt2, paaa, pbbb, pccc ) =                         & 
    79          &   paaa * pt2 + zinvtwopi * pbbb * SIN(pccc + ztwopi * pt2)   & 
    80          & - paaa * pt1 - zinvtwopi * pbbb * SIN(pccc + ztwopi * pt1) 
    8181      !!--------------------------------------------------------------------- 
    8282      ! 
    8383      ! Initialization 
    8484      ! -------------- 
    85       ztwopi    = 2._wp * rpi 
    86       zinvtwopi = 1._wp / ztwopi 
    87       zconvrad  = ztwopi / 360._wp 
    88  
    8985      ! When are we during the day (from 0 to 1) 
    9086      zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdt ) / rday 
    9187      zup = zlo + ( REAL(nn_fsbc, wp)     * rdt ) / rday 
    92       !                                           
    93       IF( nday_qsr == -1 ) THEN       ! first time step only   
     88      ! 
     89      IF( nday_qsr == -1 ) THEN       ! first time step only 
    9490         IF(lwp) THEN 
    9591            WRITE(numout,*) 
    … …  
    9894            WRITE(numout,*) 
    9995         ENDIF 
     96      END IF 
     97       
     98      ! Setting parameters for each new day: 
     99      CALL sbc_dcy_param() 
     100 
     101      !CALL iom_put( "rdusk_dcy", rdusk_dcy(:,:)*tmask(:,:,1) ) !LB 
     102      !CALL iom_put( "rdawn_dcy", rdawn_dcy(:,:)*tmask(:,:,1) ) !LB 
     103      !CALL iom_put( "rscal_dcy", rscal(:,:)*tmask(:,:,1) ) !LB 
     104 
     105 
     106      !     3. update qsr with the diurnal cycle 
     107      !     ------------------------------------ 
     108 
     109      imask_night(:,:) = 0 
     110      DO jj = 1, jpj 
     111         DO ji = 1, jpi 
     112            ztmpm = 0._wp 
     113            IF( ABS(rab(ji,jj)) < 1. ) THEN         ! day duration is less than 24h 
     114               ! 
     115               IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN       ! day time in one part 
     116                  zlousd = MAX(zlo, rdawn_dcy(ji,jj)) 
     117                  zlousd = MIN(zlousd, zup) 
     118                  zupusd = MIN(zup, rdusk_dcy(ji,jj)) 
     119                  zupusd = MAX(zupusd, zlo) 
     120                  ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     121                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
     122                  ztmpm = zupusd - zlousd 
     123                  IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 
     124                  ! 
     125               ELSE                                         ! day time in two parts 
     126                  zlousd = MIN(zlo, rdusk_dcy(ji,jj)) 
     127                  zupusd = MIN(zup, rdusk_dcy(ji,jj)) 
     128                  ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     129                  ztmpm1=zupusd-zlousd 
     130                  zlousd = MAX(zlo, rdawn_dcy(ji,jj)) 
     131                  zupusd = MAX(zup, rdawn_dcy(ji,jj)) 
     132                  ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     133                  ztmpm2 =zupusd-zlousd 
     134                  ztmp = ztmp1 + ztmp2 
     135                  ztmpm = ztmpm1 + ztmpm2 
     136                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
     137                  IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1 
     138               ENDIF 
     139            ELSE                                   ! 24h light or 24h night 
     140               ! 
     141               IF( raa(ji,jj) > rbb(ji,jj) ) THEN           ! 24h day 
     142                  ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     143                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
     144                  imask_night(ji,jj) = 0 
     145                  ! 
     146               ELSE                                         ! No day 
     147                  zqsrout(ji,jj) = 0.0_wp 
     148                  imask_night(ji,jj) = 1 
     149               ENDIF 
     150            ENDIF 
     151         END DO 
     152      END DO 
     153      ! 
     154      IF( PRESENT(l_mask) .AND. l_mask ) THEN 
     155         zqsrout(:,:) = float(imask_night(:,:)) 
     156      ENDIF 
     157      ! 
     158   END FUNCTION sbc_dcy 
     159 
     160 
     161   SUBROUTINE sbc_dcy_param( ) 
     162      !! 
     163      INTEGER  ::   ji, jj                                       ! dummy loop indices 
     164      !INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask 
     165      REAL(wp) ::   zdsws, zdecrad, ztx, zsin, zcos 
     166      REAL(wp) ::   ztmp, ztest 
     167      !---------------------------statement functions------------------------ 
     168      ! 
     169      IF( nday_qsr == -1 ) THEN       ! first time step only 
    100170         ! allocate sbcdcy arrays 
    101171         IF( sbc_dcy_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc : unable to allocate arrays' ) 
    102172         ! Compute rcc needed to compute the time integral of the diurnal cycle 
    103          rcc(:,:) = zconvrad * glamt(:,:) - rpi 
     173         rcc(:,:) = rad * glamt(:,:) - rpi 
    104174         ! time of midday 
    105175         rtmd(:,:) = 0.5_wp - glamt(:,:) / 360._wp 
    106176         rtmd(:,:) = MOD( (rtmd(:,:) + 1._wp) , 1._wp) 
    107177      ENDIF 
    108  
    109       ! If this is a new day, we have to update the dawn, dusk and scaling function   
     178       
     179      ! If this is a new day, we have to update the dawn, dusk and scaling function 
    110180      !---------------------- 
    111      
    112       !     2.1 dawn and dusk   
    113  
    114       ! nday is the number of days since the beginning of the current month  
    115       IF( nday_qsr /= nday ) THEN  
     181 
     182      !     2.1 dawn and dusk 
     183 
     184      ! nday is the number of days since the beginning of the current month 
     185      IF( nday_qsr /= nday ) THEN 
    116186         ! save the day of the year and the daily mean of qsr 
    117          nday_qsr = nday  
    118          ! number of days since the previous winter solstice (supposed to be always 21 December)          
     187         nday_qsr = nday 
     188         ! number of days since the previous winter solstice (supposed to be always 21 December) 
    119189         zdsws = REAL(11 + nday_year, wp) 
    120190         ! declination of the earths orbit 
    121          zdecrad = (-23.5_wp * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) ) 
     191         zdecrad = (-23.5_wp * rad) * COS( zdsws * 2._wp*rpi / REAL(nyear_len(1),wp) ) 
    122192         ! Compute A and B needed to compute the time integral of the diurnal cycle 
    123193 
    … …  
    125195         DO jj = 1, jpj 
    126196            DO ji = 1, jpi 
    127                ztmp = zconvrad * gphit(ji,jj) 
     197               ztmp = rad * gphit(ji,jj) 
    128198               raa(ji,jj) = SIN( ztmp ) * zsin 
    129199               rbb(ji,jj) = COS( ztmp ) * zcos 
    130             END DO   
    131          END DO   
     200            END DO 
     201         END DO 
    132202         ! Compute the time of dawn and dusk 
    133203 
    134          ! rab to test if the day time is equal to 0, less than 24h of full day         
     204         ! rab to test if the day time is equal to 0, less than 24h of full day 
    135205         rab(:,:) = -raa(:,:) / rbb(:,:) 
    136206         DO jj = 1, jpj 
    137207            DO ji = 1, jpi 
    138208               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
    139          ! When is it night? 
    140                   ztx = zinvtwopi * (ACOS(rab(ji,jj)) - rcc(ji,jj)) 
    141                   ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + ztwopi * ztx ) 
    142          ! is it dawn or dusk? 
     209                  ! When is it night? 
     210                  ztx = 1._wp/(2._wp*rpi) * (ACOS(rab(ji,jj)) - rcc(ji,jj)) 
     211                  ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + 2._wp*rpi * ztx ) 
     212                  ! is it dawn or dusk? 
    143213                  IF ( ztest > 0._wp ) THEN 
    144                      rdawn(ji,jj) = ztx 
    145                      rdusk(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn(ji,jj) ) 
     214                     rdawn_dcy(ji,jj) = ztx 
     215                     rdusk_dcy(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn_dcy(ji,jj) ) 
    146216                  ELSE 
    147                      rdusk(ji,jj) = ztx 
    148                      rdawn(ji,jj) = rtmd(ji,jj) - ( rdusk(ji,jj) - rtmd(ji,jj) ) 
     217                     rdusk_dcy(ji,jj) = ztx 
     218                     rdawn_dcy(ji,jj) = rtmd(ji,jj) - ( rdusk_dcy(ji,jj) - rtmd(ji,jj) ) 
    149219                  ENDIF 
    150220               ELSE 
    151                   rdawn(ji,jj) = rtmd(ji,jj) + 0.5_wp 
    152                   rdusk(ji,jj) = rdawn(ji,jj) 
    153                ENDIF 
    154              END DO   
    155          END DO   
    156          rdawn(:,:) = MOD( (rdawn(:,:) + 1._wp), 1._wp ) 
    157          rdusk(:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp ) 
     221                  rdawn_dcy(ji,jj) = rtmd(ji,jj) + 0.5_wp 
     222                  rdusk_dcy(ji,jj) = rdawn_dcy(ji,jj) 
     223               ENDIF 
     224            END DO 
     225         END DO 
     226         rdawn_dcy(:,:) = MOD( (rdawn_dcy(:,:) + 1._wp), 1._wp ) 
     227         rdusk_dcy(:,:) = MOD( (rdusk_dcy(:,:) + 1._wp), 1._wp ) 
    158228         !     2.2 Compute the scaling function: 
    159229         !         S* = the inverse of the time integral of the diurnal cycle from dawn to dusk 
    … …  
    164234               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
    165235                  rscal(ji,jj) = 0.0_wp 
    166                   IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN      ! day time in one part 
    167                      IF( (rdusk(ji,jj) - rdawn(ji,jj) ) .ge. 0.001_wp ) THEN 
    168                        rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    169                        rscal(ji,jj) = 1._wp / rscal(ji,jj) 
     236                  IF ( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN      ! day time in one part 
     237                     IF( (rdusk_dcy(ji,jj) - rdawn_dcy(ji,jj) ) .ge. 0.001_wp ) THEN 
     238                        rscal(ji,jj) = fintegral(rdawn_dcy(ji,jj), rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     239                        rscal(ji,jj) = 1._wp / rscal(ji,jj) 
    170240                     ENDIF 
    171241                  ELSE                                         ! day time in two parts 
    172                      IF( (rdusk(ji,jj) + (1._wp - rdawn(ji,jj)) ) .ge. 0.001_wp ) THEN 
    173                        rscal(ji,jj) = fintegral(0._wp, rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   & 
    174                           &         + fintegral(rdawn(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    175                        rscal(ji,jj) = 1. / rscal(ji,jj) 
     242                     IF( (rdusk_dcy(ji,jj) + (1._wp - rdawn_dcy(ji,jj)) ) .ge. 0.001_wp ) THEN 
     243                        rscal(ji,jj) = fintegral(0._wp, rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   & 
     244                           &         + fintegral(rdawn_dcy(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     245                        rscal(ji,jj) = 1. / rscal(ji,jj) 
    176246                     ENDIF 
    177247                  ENDIF 
    178248               ELSE 
    179249                  IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day 
    180                      rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
     250                     rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
    181251                     rscal(ji,jj) = 1._wp / rscal(ji,jj) 
    182252                  ELSE                                          ! No day 
    … …  
    184254                  ENDIF 
    185255               ENDIF 
    186             END DO   
    187          END DO   
     256            END DO 
     257         END DO 
    188258         ! 
    189259         ztmp = rday / ( rdt * REAL(nn_fsbc, wp) ) 
    190260         rscal(:,:) = rscal(:,:) * ztmp 
    191261         ! 
    192       ENDIF  
    193          !     3. update qsr with the diurnal cycle 
    194          !     ------------------------------------ 
    195  
    196       imask_night(:,:) = 0 
    197       DO jj = 1, jpj 
    198          DO ji = 1, jpi 
    199             ztmpm = 0._wp 
    200             IF( ABS(rab(ji,jj)) < 1. ) THEN         ! day duration is less than 24h 
    201                ! 
    202                IF( rdawn(ji,jj) < rdusk(ji,jj) ) THEN       ! day time in one part 
    203                   zlousd = MAX(zlo, rdawn(ji,jj)) 
    204                   zlousd = MIN(zlousd, zup) 
    205                   zupusd = MIN(zup, rdusk(ji,jj)) 
    206                   zupusd = MAX(zupusd, zlo) 
    207                   ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    208                   zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
    209                   ztmpm = zupusd - zlousd 
    210                   IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 
    211                   ! 
    212                ELSE                                         ! day time in two parts 
    213                   zlousd = MIN(zlo, rdusk(ji,jj)) 
    214                   zupusd = MIN(zup, rdusk(ji,jj)) 
    215                   ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    216                   ztmpm1=zupusd-zlousd 
    217                   zlousd = MAX(zlo, rdawn(ji,jj)) 
    218                   zupusd = MAX(zup, rdawn(ji,jj)) 
    219                   ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    220                   ztmpm2 =zupusd-zlousd 
    221                   ztmp = ztmp1 + ztmp2 
    222                   ztmpm = ztmpm1 + ztmpm2 
    223                   zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
    224                   IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1 
    225                ENDIF 
    226             ELSE                                   ! 24h light or 24h night 
    227                ! 
    228                IF( raa(ji,jj) > rbb(ji,jj) ) THEN           ! 24h day 
    229                   ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    230                   zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
    231                   imask_night(ji,jj) = 0 
    232                   ! 
    233                ELSE                                         ! No day 
    234                   zqsrout(ji,jj) = 0.0_wp 
    235                   imask_night(ji,jj) = 1 
    236                ENDIF 
    237             ENDIF 
    238          END DO   
    239       END DO   
    240       ! 
    241       IF( PRESENT(l_mask) .AND. l_mask ) THEN 
    242          zqsrout(:,:) = float(imask_night(:,:)) 
    243       ENDIF 
    244       ! 
    245    END FUNCTION sbc_dcy 
     262      ENDIF !IF( nday_qsr /= nday ) 
     263      ! 
     264   END SUBROUTINE sbc_dcy_param 
     265 
     266 
     267   FUNCTION fintegral( pt1, pt2, paaa, pbbb, pccc ) 
     268      REAL(wp), INTENT(in) :: pt1, pt2, paaa, pbbb, pccc 
     269      REAL(wp) :: fintegral 
     270      fintegral =   paaa * pt2 + 1._wp/(2._wp*rpi) * pbbb * SIN(pccc + 2._wp*rpi*pt2)   & 
     271         &        - paaa * pt1 - 1._wp/(2._wp*rpi) * pbbb * SIN(pccc + 2._wp*rpi*pt1) 
     272   END FUNCTION fintegral 
    246273 
    247274   !!====================================================================== 

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

    r10499 r11672  
    254254 
    255255      !                          ! Choice of the Surface Boudary Condition (set nsbc) 
     256      nday_qsr = -1   ! allow initialization at the 1st call !LB: now warm-layer of COARE* calls "sbc_dcy_param" of sbcdcy.F90! 
    256257      IF( ln_dm2dc ) THEN           !* daily mean to diurnal cycle 
    257          nday_qsr = -1   ! allow initialization at the 1st call 
     258         !LB:nday_qsr = -1   ! allow initialization at the 1st call 
    258259         IF( .NOT.( ln_flx .OR. ln_blk ) .AND. nn_components /= jp_iam_opa )   & 
    259260            &   CALL ctl_stop( 'qsr diurnal cycle from daily values requires a flux or bulk formulation' ) 
Note: See TracChangeset for help on using the changeset viewer.