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

Ignore:
Timestamp:
2019-08-16T12:32:43+02:00 (5 years ago)
Author:
mattmartin
Message:

Introduction of stochastic physics in NEMO, based on Andrea Storto's code.
For details, see ticket https://code.metoffice.gov.uk/trac/utils/ticket/251.

File:
1 edited

Legend:

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

    r9288 r11442  
    5151   USE par_ice_2 
    5252#endif 
     53   USE stopack 
    5354 
    5455   IMPLICIT NONE 
     
    8990   REAL(wp) ::   rn_efac     ! multiplication factor for evaporation (clem) 
    9091   REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 
     92   REAL(wp), ALLOCATABLE, SAVE ::   rn_vfac0(:,:) ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 
    9193   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements 
    9294   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements 
     
    198200         sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 
    199201         ! 
     202         ALLOCATE ( rn_vfac0(jpi,jpj) ) 
     203         rn_vfac0(:,:) = rn_vfac 
     204         ! 
     205      ENDIF 
     206 
     207      IF( ln_stopack .AND. nn_spp_relw > 0 ) THEN 
     208         rn_vfac0(:,:) = rn_vfac 
     209         CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw ) 
    200210      ENDIF 
    201211 
     
    291301      DO jj = 2, jpjm1 
    292302         DO ji = fs_2, fs_jpim1   ! vect. opt. 
    293             zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
    294             zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
     303            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
     304            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
    295305         END DO 
    296306      END DO 
     
    467477            DO ji = 2, jpim1   ! B grid : NO vector opt 
    468478               ! ... scalar wind at I-point (fld being at T-point) 
    469                zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
    470                   &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
    471                zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    472                   &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
     479               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)    & 
     480                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) & 
     481                  &      - rn_vfac0(ji,jj) * u_ice(ji,jj) 
     482               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)    & 
     483                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) & 
     484                  &      - rn_vfac0(ji,jj) * v_ice(ji,jj) 
    473485               zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    474486               ! ... ice stress at I-point 
     
    476488               vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
    477489               ! ... scalar wind at T-point (fld being at T-point) 
    478                zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
    479                   &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  ) 
    480                zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
    481                   &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
     490               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1)                                         & 
     491                  &      - rn_vfac0(ji,jj) * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
     492                  &                                  + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  ) 
     493               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1)                                         & 
     494                  &      - rn_vfac0(ji,jj) * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
     495                  &                                  + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
    482496               wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    483497            END DO 
     
    490504         DO jj = 2, jpj 
    491505            DO ji = fs_2, jpi   ! vect. opt. 
    492                zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
    493                zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
     506               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
     507               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
    494508               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    495509            END DO 
     
    497511         DO jj = 2, jpjm1 
    498512            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    499                utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
    500                   &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    501                vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
    502                   &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
     513               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )             & 
     514                  &              * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) & 
     515                  &                  - rn_vfac0(ji,jj) * u_ice(ji,jj) ) 
     516               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )             & 
     517                  &              * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) & 
     518                  &                  - rn_vfac0(ji,jj) * v_ice(ji,jj) ) 
    503519            END DO 
    504520         END DO 
     
    645661      DO jl = 1, jpl 
    646662         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 
    647                                    ! But we do not have Tice => consider it at 0°C => evap=0  
     663                                   ! But we do not have Tice => consider it at 0 degC => evap=0  
    648664      END DO 
    649665 
Note: See TracChangeset for help on using the changeset viewer.