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

Changeset 11217


Ignore:
Timestamp:
2019-07-04T17:31:44+02:00 (5 years ago)
Author:
laurent
Message:

LB: more readable "cool-skin/warm-layer" organisation + no "cool-skin/warm-layer" scheme in the presence of sea-ice.

File:
1 edited

Legend:

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

    r11215 r11217  
    1717   !!                 !                        ==> based on AeroBulk (https://github.com/brodeau/aerobulk/) 
    1818   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine 
    19    !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle)  
     19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle) 
    2020   !!---------------------------------------------------------------------- 
    2121 
     
    2424   !!   sbc_blk       : bulk formulation as ocean surface boundary condition 
    2525   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean 
    26    !!             sea-ice case only :  
     26   !!             sea-ice case only : 
    2727   !!   blk_ice_tau   : provide the air-ice stress 
    2828   !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface 
    2929   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 
    3030   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
    31    !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
     31   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
    3232   !!---------------------------------------------------------------------- 
    3333   USE oce            ! ocean dynamics and tracers 
     
    4545   USE icethd_dh      ! for CALL ice_thd_snwblow 
    4646#endif 
    47    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
    48    USE sbcblk_algo_coare    ! => turb_coare    : COAREv3.0 (Fairall et al. 2003)  
     47   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
     48   USE sbcblk_algo_coare    ! => turb_coare    : COAREv3.0 (Fairall et al. 2003) 
    4949   USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013) 
    50    USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31)  
     50   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31) 
    5151   ! 
    5252   USE iom            ! I/O manager library 
     
    5858   USE sbcblk_phy     !LB: all thermodynamics functions in the marine boundary layer, rho_air, q_sat, etc... 
    5959 
    60     
     60 
    6161   IMPLICIT NONE 
    6262   PRIVATE 
     
    6868   PUBLIC   blk_ice_flx   ! routine called in icesbc 
    6969   PUBLIC   blk_ice_qcn   ! routine called in icesbc 
    70 #endif  
     70#endif 
    7171 
    7272   INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read 
     
    9696   REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements 
    9797   REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements 
    98 !!gm ref namelist initialize it so remove the setting to false below 
     98   !!gm ref namelist initialize it so remove the setting to false below 
    9999   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012) 
    100100   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) 
     
    111111   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature if .true. (specific humidity espected if .false.) !LB 
    112112   !LB. 
    113     
     113 
    114114   INTEGER  ::   nblk           ! choice of the bulk algorithm 
    115115   !                            ! associated indices: 
     
    159159      !! ** Purpose :   choose and initialize a bulk formulae formulation 
    160160      !! 
    161       !! ** Method  :  
     161      !! ** Method  : 
    162162      !! 
    163163      !!---------------------------------------------------------------------- 
     
    173173         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                & 
    174174         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
    175          &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
     175         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
    176176         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           & 
    177177         &                 ln_skin, ln_humi_dpt                                                        ! cool-skin / warm-layer param. !LB 
     
    181181      IF( sbc_blk_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 
    182182      ! 
    183       !                             !** read bulk namelist   
     183      !                             !** read bulk namelist 
    184184      REWIND( numnam_ref )                !* Namelist namsbc_blk in reference namelist : bulk parameters 
    185185      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) 
     
    194194      !                             !** initialization of the chosen bulk formulae (+ check) 
    195195      !                                   !* select the bulk chosen in the namelist and check the choice 
    196                                                                ioptio = 0 
    197       IF( ln_NCAR      ) THEN   ;   nblk =  np_NCAR        ;   ioptio = ioptio + 1   ;   ENDIF 
    198       IF( ln_COARE_3p0 ) THEN   ;   nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1   ;   ENDIF 
    199       IF( ln_COARE_3p5 ) THEN   ;   nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1   ;   ENDIF 
    200       IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF 
     196      ioptio = 0 
     197      IF( ln_NCAR      ) THEN 
     198         nblk =  np_NCAR        ;   ioptio = ioptio + 1 
     199      ENDIF 
     200      IF( ln_COARE_3p0 ) THEN 
     201         nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1 
     202      ENDIF 
     203      IF( ln_COARE_3p5 ) THEN 
     204         nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1 
     205      ENDIF 
     206      IF( ln_ECMWF     ) THEN 
     207         nblk =  np_ECMWF       ;   ioptio = ioptio + 1 
     208      ENDIF 
    201209      ! 
    202210      !LB: 
     
    208216      END IF 
    209217      !LB. 
    210           
     218 
    211219      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 
    212220      ! 
    213221      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr 
    214222         IF( sn_qsr%nfreqh /= 24 )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' ) 
    215          IF( sn_qsr%ln_tint ) THEN  
     223         IF( sn_qsr%ln_tint ) THEN 
    216224            CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   & 
    217225               &           '              ==> We force time interpolation = .false. for qsr' ) 
    218             sn_qsr%ln_tint = .false. 
     226            sn_qsr%ln_tint = .FALSE. 
    219227         ENDIF 
    220228      ENDIF 
     
    245253      ! 
    246254      IF ( ln_wave ) THEN 
    247       !Activated wave module but neither drag nor stokes drift activated 
     255         !Activated wave module but neither drag nor stokes drift activated 
    248256         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN 
    249257            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' ) 
    250       !drag coefficient read from wave model definable only with mfs bulk formulae and core  
    251          ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN        
    252              CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 
     258            !drag coefficient read from wave model definable only with mfs bulk formulae and core 
     259         ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN 
     260            CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 
    253261         ELSEIF (ln_stcor .AND. .NOT. ln_sdw)                             THEN 
    254              CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
     262            CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
    255263         ENDIF 
    256264      ELSE 
    257       IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                &  
    258          &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
    259          &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
    260          &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
    261          &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      &   
    262          &                  'or Stokes-Coriolis term (ln_stcori=T)'  ) 
    263       ENDIF  
    264       ! 
    265       !            
     265         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                & 
     266            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
     267            &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
     268            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
     269            &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      & 
     270            &                  'or Stokes-Coriolis term (ln_stcori=T)'  ) 
     271      ENDIF 
     272      ! 
     273      ! 
    266274      IF(lwp) THEN                     !** Control print 
    267275         ! 
    268          WRITE(numout,*)                  !* namelist  
     276         WRITE(numout,*)                  !* namelist 
    269277         WRITE(numout,*) '   Namelist namsbc_blk (other than data information):' 
    270278         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR 
     
    344352      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN 
    345353         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1) 
    346          IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 
    347          ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)  
    348          ENDIF  
     354         IF( ln_dm2dc ) THEN 
     355            qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 
     356         ELSE 
     357            qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
     358         ENDIF 
    349359         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1) 
    350360         !LB: 
     
    411421 
    412422      ! ... components ( U10m - U_oce ) at T-point (unmasked) 
    413 !!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ??? 
     423      !!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ??? 
    414424      zwnd_i(:,:) = 0._wp 
    415425      zwnd_j(:,:) = 0._wp 
     
    440450      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave 
    441451      zztmp = 1. - albo 
    442       IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 
    443       ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
     452      IF( ln_dm2dc ) THEN 
     453         qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 
     454      ELSE 
     455         qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
    444456      ENDIF 
    445457 
     
    451463      ! ... specific humidity at SST and IST tmask( 
    452464      zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 
     465 
     466      IF ( ln_skin ) THEN 
     467         zqsb(:,:) = zst(:,:) !LB: using array zqsb to backup original zst before skin action 
     468         zqla(:,:) = zsq(:,:) !LB: using array zqla to backup original zsq before skin action 
     469      END IF 
    453470 
    454471      !LB: 
     
    458475         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of d_air and slp !' 
    459476         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    460       ELSE          
     477      ELSE 
    461478         zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity! 
    462479      END IF 
     
    467484      !!    (since reanalysis products provide T at z, not theta !) 
    468485      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt 
    469  
    470       SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    471       ! 
    472       CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! NCAR-COREv2 
    473          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    474  
    475       !LB: Skin!!!          
    476       CASE( np_COARE_3p0 ) 
    477          IF ( ln_skin ) THEN             
     486       
     487 
     488      IF ( ln_skin ) THEN 
     489          
     490         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
     491 
     492         CASE( np_COARE_3p0 ) 
    478493            CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0 
    479494               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 
    480495               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1)) 
     496 
     497         CASE( np_ECMWF     ) 
     498            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF 
     499               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 
     500               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1)) 
     501 
     502         CASE DEFAULT 
     503            CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin==.true."' ) 
     504         END SELECT 
     505 
     506         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and zsq: 
     507         WHERE ( fr_i < 0.001_wp ) 
     508            ! zst and zsq have been updated by cool-skin/warm-layer scheme and we keep it!!! 
    481509            zst(:,:) = zst(:,:)*tmask(:,:,1) 
    482510            zsq(:,:) = zsq(:,:)*tmask(:,:,1) 
    483          ELSE 
    484             IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare" WITHOUT CSWL optional arrays!!!'             
     511         ELSEWHERE 
     512            ! we forget about the update... 
     513            zst(:,:) = zqsb(:,:) !LB: using what we backed up before skin-algo 
     514            zsq(:,:) = zqla(:,:) !LB:  "   "   " 
     515         END WHERE 
     516 
     517         !LB: Update of tsk, the official array for skin temperature 
     518         tsk(:,:) = zst(:,:) 
     519 
     520      ELSE !IF ( ln_skin ) 
     521 
     522          
     523         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
     524            ! 
     525         CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! NCAR-COREv2 
     526            &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     527 
     528         CASE( np_COARE_3p0 ) 
     529            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare" WITHOUT CSWL optional arrays!!!' 
    485530            CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0 
    486531               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    487          END IF 
    488              
    489       CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.5 
    490          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    491  
    492       !LB: Skin!!! 
    493       CASE( np_ECMWF     ) 
    494          IF ( ln_skin ) THEN             
    495             CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF          
    496                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 
    497                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1)) 
    498             zst(:,:) = zst(:,:)*tmask(:,:,1) 
    499             zsq(:,:) = zsq(:,:)*tmask(:,:,1) 
    500          ELSE 
     532 
     533         CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.5 
     534            &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     535 
     536         CASE( np_ECMWF     ) 
    501537            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!' 
    502             CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF          
     538            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF 
    503539               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    504          END IF 
    505          !LB. 
    506           
    507       CASE DEFAULT 
    508          CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
    509       END SELECT 
    510  
    511  
    512       !LB: cool-skin/warm-layer: 
    513       IF ( ln_skin ) tsk(:,:) = zst(:,:) 
     540 
     541         CASE DEFAULT 
     542            CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
     543         END SELECT 
     544 
     545      END IF ! IF ( ln_skin ) 
     546 
     547 
    514548 
    515549 
     
    521555      END IF 
    522556 
    523 !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
    524 !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef. 
     557      !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
     558      !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef. 
    525559 
    526560      DO jj = 1, jpj             ! tau module, i and j component 
     
    612646      ! 
    613647      !!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN 
    614          CALL iom_put( "rho_air" ,   rhoa )                 ! output air density (kg/m^3) !#LB 
    615          CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean 
    616          CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean 
    617          CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean 
    618          CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean 
    619          CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean 
    620          CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean 
    621          CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean 
    622          tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 
    623          sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 
    624          CALL iom_put( 'snowpre', sprecip )                 ! Snow 
    625          CALL iom_put( 'precip' , tprecip )                 ! Total precipitation 
    626          IF ( ln_skin ) THEN 
    627             CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius 
    628             CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference... 
    629          END IF 
     648      CALL iom_put( "rho_air" ,   rhoa )                 ! output air density (kg/m^3) !#LB 
     649      CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean 
     650      CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean 
     651      CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean 
     652      CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean 
     653      CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean 
     654      CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean 
     655      CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean 
     656      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 
     657      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 
     658      CALL iom_put( 'snowpre', sprecip )                 ! Snow 
     659      CALL iom_put( 'precip' , tprecip )                 ! Total precipitation 
     660      IF ( ln_skin ) THEN 
     661         CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius 
     662         CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference... 
     663      END IF 
    630664      !!#LB. ENDIF 
    631665      ! 
     
    650684   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 
    651685   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
    652    !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
     686   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
    653687   !!---------------------------------------------------------------------- 
    654688 
     
    693727         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical 
    694728      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations 
    695          CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )  
    696       ENDIF 
    697  
    698 !!      CALL iom_put( "rCd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef. 
    699 !!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef. 
     729         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 
     730      ENDIF 
     731 
     732      !!      CALL iom_put( "rCd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef. 
     733      !!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef. 
    700734 
    701735      ! local scalars ( place there for vector optimisation purposes) 
     
    776810 
    777811      zztmp = 1. / ( 1. - albo ) 
    778       WHERE( ptsu(:,:,:) /= 0._wp )   ;   z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 
    779       ELSEWHERE                       ;   z1_st(:,:,:) = 0._wp 
     812      WHERE( ptsu(:,:,:) /= 0._wp ) 
     813         z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 
     814      ELSEWHERE 
     815         z1_st(:,:,:) = 0._wp 
    780816      END WHERE 
    781817      !                                     ! ========================== ! 
     
    866902      DO jl = 1, jpl 
    867903         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 
    868          !                         ! But we do not have Tice => consider it at 0degC => evap=0  
     904         !                         ! But we do not have Tice => consider it at 0degC => evap=0 
    869905      END DO 
    870906 
     
    873909      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1 
    874910      ! 
    875       WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
     911      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm 
    876912         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
    877913      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm 
    878914         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 
    879915      ELSEWHERE                                                         ! zero when hs>0 
    880          qtr_ice_top(:,:,:) = 0._wp  
     916         qtr_ice_top(:,:,:) = 0._wp 
    881917      END WHERE 
    882918      ! 
     
    891927      ! 
    892928   END SUBROUTINE blk_ice_flx 
    893     
     929 
    894930 
    895931   SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi ) 
     
    900936      !!                to force sea ice / snow thermodynamics 
    901937      !!                in the case conduction flux is emulated 
    902       !!                 
     938      !! 
    903939      !! ** Method  :   compute surface energy balance assuming neglecting heat storage 
    904940      !!                following the 0-layer Semtner (1976) approach 
     
    925961      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor 
    926962      !!--------------------------------------------------------------------- 
    927        
     963 
    928964      ! -------------------------------------! 
    929965      !      I   Enhanced conduction factor  ! 
     
    933969      ! 
    934970      zgfac(:,:,:) = 1._wp 
    935        
     971 
    936972      IF( ld_virtual_itd ) THEN 
    937973         ! 
     
    939975         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 
    940976         zfac3 = 2._wp / zepsilon 
    941          !    
    942          DO jl = 1, jpl                 
     977         ! 
     978         DO jl = 1, jpl 
    943979            DO jj = 1 , jpj 
    944980               DO ji = 1, jpi 
     
    948984            END DO 
    949985         END DO 
    950          !       
    951       ENDIF 
    952        
     986         ! 
     987      ENDIF 
     988 
    953989      ! -------------------------------------------------------------! 
    954990      !      II   Surface temperature and conduction flux            ! 
     
    960996         DO jj = 1 , jpj 
    961997            DO ji = 1, jpi 
    962                !                     
     998               ! 
    963999               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness 
    9641000                  &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 
     
    9771013               qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 
    9781014               qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) )  & 
    979                              &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 
     1015                  &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 
    9801016 
    9811017               ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 
    982                hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl)  
     1018               hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 
    9831019 
    9841020            END DO 
    9851021         END DO 
    9861022         ! 
    987       END DO  
    988       !       
     1023      END DO 
     1024      ! 
    9891025   END SUBROUTINE blk_ice_qcn 
    990     
     1026 
    9911027 
    9921028   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
     
    9941030      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
    9951031      !! 
    996       !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m  
     1032      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m 
    9971033      !!                 to make it dependent on edges at leads, melt ponds and flows. 
    9981034      !!                 After some approximations, this can be resumed to a dependency 
    9991035      !!                 on ice concentration. 
    1000       !!                 
     1036      !! 
    10011037      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50) 
    10021038      !!                 with the highest level of approximation: level4, eq.(59) 
     
    10101046      !! 
    10111047      !!                 This new drag has a parabolic shape (as a function of A) starting at 
    1012       !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5  
     1048      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 
    10131049      !!                 and going down to Cdi(say 1.4e-3) for A=1 
    10141050      !! 
     
    10371073      Cd(:,:) = rCd_ice +  &                                                         ! pure ice drag 
    10381074         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology 
    1039        
     1075 
    10401076   END SUBROUTINE Cdn10_Lupkes2012 
    10411077 
     
    10461082      !! 
    10471083      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation 
    1048       !!                 between sea-ice and atmosphere with distinct momentum  
    1049       !!                 and heat coefficients depending on sea-ice concentration  
     1084      !!                 between sea-ice and atmosphere with distinct momentum 
     1085      !!                 and heat coefficients depending on sea-ice concentration 
    10501086      !!                 and atmospheric stability (no meltponds effect for now). 
    1051       !!                 
     1087      !! 
    10521088      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015) 
    10531089      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 
    10541090      !!                 it considers specific skin and form drags (Andreas et al. 2010) 
    1055       !!                 to compute neutral transfert coefficients for both heat and  
     1091      !!                 to compute neutral transfert coefficients for both heat and 
    10561092      !!                 momemtum fluxes. Atmospheric stability effect on transfert 
    10571093      !!                 coefficient is also taken into account following Louis (1979). 
     
    10931129 
    10941130      ! mean temperature 
    1095       WHERE( at_i_b(:,:) > 1.e-20 )   ;   ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:) 
    1096       ELSEWHERE                       ;   ztm_su(:,:) = rt0 
     1131      WHERE( at_i_b(:,:) > 1.e-20 ) 
     1132         ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:) 
     1133      ELSEWHERE 
     1134         ztm_su(:,:) = rt0 
    10971135      ENDWHERE 
    1098        
     1136 
    10991137      ! Momentum Neutral Transfert Coefficients (should be a constant) 
    11001138      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40 
     
    11051143      ! Heat Neutral Transfert Coefficients 
    11061144      zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) )   ! Eq. 50 + Eq. 52 (cf Lupkes email for details) 
    1107       
     1145 
    11081146      ! Atmospheric and Surface Variables 
    11091147      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin 
     
    11171155            zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice 
    11181156            zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu 
    1119              
     1157 
    11201158            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 
    11211159            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean 
    11221160            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice 
    1123              
     1161 
    11241162            ! Momentum and Heat Neutral Transfert Coefficients 
    11251163            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40 
    1126             zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53  
    1127                         
     1164            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53 
     1165 
    11281166            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 
    11291167            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 
     
    11371175               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28 
    11381176            ENDIF 
    1139              
     1177 
    11401178            IF( zrib_i <= 0._wp ) THEN 
    11411179               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9 
     
    11451183               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27 
    11461184            ENDIF 
    1147              
     1185 
    11481186            ! Momentum Transfert Coefficients (Eq. 38) 
    11491187            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  & 
    11501188               &        zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 
    1151              
     1189 
    11521190            ! Heat Transfert Coefficients (Eq. 49) 
    11531191            Ch(ji,jj) = zChn_skin_ice *   zfhi +  & 
Note: See TracChangeset for help on using the changeset viewer.