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 13088 for branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 – NEMO

Ignore:
Timestamp:
2020-06-10T13:13:39+02:00 (4 years ago)
Author:
jwhile
Message:

Bug fixes for 1D running

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r11442 r13088  
    161161         !                                      ! ====================== ! 
    162162         ! 
    163          rn_sfac = 1._wp       ! Default to one if missing from namelist  
     163         rn_sfac = 1._wp       ! Default to one if missing from namelist 
    164164         REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters 
    165165         READ  ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901) 
     
    205205      ENDIF 
    206206 
     207#if defined key_traldf_c2d || key_traldf_c3d 
    207208      IF( ln_stopack .AND. nn_spp_relw > 0 ) THEN 
    208209         rn_vfac0(:,:) = rn_vfac 
    209210         CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw ) 
    210211      ENDIF 
     212#else 
     213      IF ( ln_stopack .AND. nn_spp_relw > 0 ) & 
     214         & CALL ctl_stop( 'sbc_blk_core: parameter perturbation will only work with '// & 
     215                          'key_traldf_c2d or key_traldf_c3d') 
     216#endif 
    211217 
    212218      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step 
     
    217223#if defined key_cice 
    218224      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN 
    219          qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1)  
     225         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1) 
    220226         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 
    221227         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
    222228         ENDIF 
    223          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)          
     229         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1) 
    224230         qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1) 
    225231         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
     
    231237      ! 
    232238   END SUBROUTINE sbc_blk_core 
    233     
    234     
     239 
     240 
    235241   SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv ) 
    236242      !!--------------------------------------------------------------------- 
     
    242248      !! ** Method  :   CORE bulk formulea for the ocean using atmospheric 
    243249      !!      fields read in sbc_read 
    244       !!  
     250      !! 
    245251      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2) 
    246252      !!              - vtau    : j-component of the stress at V-point  (N/m2) 
     
    280286      ! local scalars ( place there for vector optimisation purposes) 
    281287      zcoef_qsatw = 0.98 * 640380. / rhoa 
    282        
     288 
    283289      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    284290 
     
    288294 
    289295      ! ... components ( U10m - U_oce ) at T-point (unmasked) 
    290       zwnd_i(:,:) = 0.e0   
     296      zwnd_i(:,:) = 0.e0 
    291297      zwnd_j(:,:) = 0.e0 
    292298#if defined key_cyclone 
     
    332338      CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   & 
    333339         &               Cd, Ch, Ce, zt_zu, zq_zu ) 
    334      
     340 
    335341      ! ... tau module, i and j component 
    336342      DO jj = 1, jpj 
     
    363369      CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
    364370 
    365      
     371 
    366372      !  Turbulent fluxes over ocean 
    367373      ! ----------------------------- 
     
    388394         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce_core: zst    : ') 
    389395      ENDIF 
    390         
     396 
    391397      ! ----------------------------------------------------------------------------- ! 
    392398      !     III    Total FLUXES                                                       ! 
     
    396402         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    397403      ! 
    398       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar  
     404      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
    399405         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip 
    400406         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
     
    437443      ! 
    438444   END SUBROUTINE blk_oce_core 
    439   
    440     
     445 
     446 
    441447#if defined key_lim2 || defined key_lim3 
    442448   SUBROUTINE blk_ice_core_tau 
     
    531537 
    532538      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau') 
    533        
     539 
    534540   END SUBROUTINE blk_ice_core_tau 
    535541 
     
    544550      !!                between atmosphere and sea-ice using CORE bulk 
    545551      !!                formulea, ice variables and read atmmospheric fields. 
    546       !!  
     552      !! 
    547553      !! caution : the net upward water flux has with mm/day unit 
    548554      !!--------------------------------------------------------------------- 
     
    564570      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_flx') 
    565571      ! 
    566       CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )  
     572      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    567573 
    568574      ! local scalars ( place there for vector optimisation purposes) 
     
    587593               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    588594               ! lw sensitivity 
    589                z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                                
     595               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 
    590596 
    591597               ! ----------------------------! 
     
    597603               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
    598604               ! Latent Heat 
    599                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                            
     605               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   & 
    600606                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    601607              ! Latent heat sensitivity for ice (Dqla/Dt) 
     
    628634 
    629635#if defined  key_lim3 
    630       CALL wrk_alloc( jpi,jpj, zevap, zsnw )  
     636      CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 
    631637 
    632638      ! --- evaporation --- ! 
     
    638644      ! --- evaporation minus precipitation --- ! 
    639645      zsnw(:,:) = 0._wp 
    640       CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing  
     646      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing 
    641647      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
    642648      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
     
    661667      DO jl = 1, jpl 
    662668         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 
    663                                    ! But we do not have Tice => consider it at 0 degC => evap=0  
     669                                   ! But we do not have Tice => consider it at 0 degC => evap=0 
    664670      END DO 
    665671 
    666       CALL wrk_dealloc( jpi,jpj, zevap, zsnw )  
     672      CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
    667673#endif 
    668674 
     
    688694      ! 
    689695      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx') 
    690        
     696 
    691697   END SUBROUTINE blk_ice_core_flx 
    692698#endif 
     
    701707      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
    702708      !! 
    703       !! ** Method : Monin Obukhov Similarity Theory  
     709      !! ** Method : Monin Obukhov Similarity Theory 
    704710      !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 
    705711      !! 
     
    711717      !!    - better first guess of stability by checking air-sea difference of virtual temperature 
    712718      !!       rather than temperature difference only... 
    713       !!    - added function "cd_neutral_10m" that uses the improved parametrization of  
     719      !!    - added function "cd_neutral_10m" that uses the improved parametrization of 
    714720      !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 
    715721      !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them 
     
    746752 
    747753      IF( nn_timing == 1 )  CALL timing_start('turb_core_2z') 
    748      
     754 
    749755      CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 
    750756      CALL wrk_alloc( jpi,jpj, zeta_u, stab ) 
     
    758764      U_zu = MAX( 0.5 , dU )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
    759765 
    760       !! First guess of stability:  
     766      !! First guess of stability: 
    761767      ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt 
    762768      stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
     
    772778      Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 
    773779      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 
    774      
     780 
    775781      !! Initializing transf. coeff. with their first guess neutral equivalents : 
    776782      Cd = ztmp0   ;   Ce = Ce_n10   ;   Ch = Ch_n10   ;   sqrt_Cd = sqrt_Cd_n10 
     
    783789         ! 
    784790         ztmp1 = T_zu - sst   ! Updating air/sea differences 
    785          ztmp2 = q_zu - q_sat  
     791         ztmp2 = q_zu - q_sat 
    786792 
    787793         ! Updating turbulent scales :   (L&Y 2004 eq. (7)) 
    788794         ztmp1  = Ch/sqrt_Cd*ztmp1    ! theta* 
    789795         ztmp2  = Ce/sqrt_Cd*ztmp2    ! q* 
    790         
     796 
    791797         ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu 
    792798 
    793799         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 
    794          ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu)  
     800         ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 
    795801         !                                                                     ( Cd*U_zu*U_zu is U*^2 at zu) 
    796802 
     
    799805         zpsi_h_u = psi_h( zeta_u ) 
    800806         zpsi_m_u = psi_m( zeta_u ) 
    801         
     807 
    802808         !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c)) 
    803809         IF ( .NOT. l_zt_equal_zu ) THEN 
     
    808814            q_zu = max(0., q_zu) 
    809815         END IF 
    810         
     816 
    811817         IF( ln_cdgw ) THEN      ! surface wave case 
    812             sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u )  
     818            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 
    813819            Cd      = sqrt_Cd * sqrt_Cd 
    814820         ELSE 
     
    820826           ztmp0 = cd_neutral_10m(ztmp0)                                                 ! Cd_n10 
    821827           sqrt_Cd_n10 = sqrt(ztmp0) 
    822         
     828 
    823829           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                     ! L&Y 2004 eq. (6b) 
    824830           stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability 
     
    827833           !! Update of transfer coefficients: 
    828834           ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)   ! L&Y 2004 eq. (10a) 
    829            Cd      = ztmp0 / ( ztmp1*ztmp1 )    
     835           Cd      = ztmp0 / ( ztmp1*ztmp1 ) 
    830836           sqrt_Cd = SQRT( Cd ) 
    831837         ENDIF 
     
    833839         ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    834840         ztmp2 = sqrt_Cd / sqrt_Cd_n10 
    835          ztmp1 = 1. + Ch_n10*ztmp0                
     841         ztmp1 = 1. + Ch_n10*ztmp0 
    836842         Ch  = Ch_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10b) 
    837843         ! 
    838          ztmp1 = 1. + Ce_n10*ztmp0                
     844         ztmp1 = 1. + Ce_n10*ztmp0 
    839845         Ce  = Ce_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    840846         ! 
     
    854860   FUNCTION cd_neutral_10m( zw10 ) 
    855861      !!---------------------------------------------------------------------- 
    856       !! Estimate of the neutral drag coefficient at 10m as a function  
     862      !! Estimate of the neutral drag coefficient at 10m as a function 
    857863      !! of neutral wind  speed at 10m 
    858864      !! 
     
    860866      !! 
    861867      !! Author: L. Brodeau, june 2014 
    862       !!----------------------------------------------------------------------     
     868      !!---------------------------------------------------------------------- 
    863869      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zw10           ! scalar wind speed at 10m (m/s) 
    864870      REAL(wp), DIMENSION(jpi,jpj)             ::   cd_neutral_10m 
    865871      ! 
    866872      REAL(wp), DIMENSION(:,:), POINTER ::   rgt33 
    867       !!----------------------------------------------------------------------     
     873      !!---------------------------------------------------------------------- 
    868874      ! 
    869875      CALL wrk_alloc( jpi,jpj, rgt33 ) 
    870876      ! 
    871877      !! When wind speed > 33 m/s => Cyclone conditions => special treatment 
    872       rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1   
     878      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1 
    873879      cd_neutral_10m = 1.e-3 * ( & 
    874880         &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. 
Note: See TracChangeset for help on using the changeset viewer.