Changeset 6389 for branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO

Timestamp:

2016-03-15T11:31:28+01:00 (8 years ago)

Author:

rfurner

Message:

Added surge flux formulation, and amended bottom friction coefficient

Location:

branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC

Files:

• sbc_oce.F90 (modified) (2 diffs)
• sbcblk_core.F90 (modified) (9 diffs)
• sbcmod.F90 (modified) (8 diffs)
• sbcsurge.F90 (added)

branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

@@ -35 +35 @@
    LOGICAL , PUBLIC ::   ln_blk_core    !: CORE bulk formulation
    LOGICAL , PUBLIC ::   ln_blk_mfs     !: MFS  bulk formulation
+   LOGICAL , PUBLIC ::   ln_sbcsurge    !: Surge model formulation
 #if defined key_oasis3
    LOGICAL , PUBLIC ::   lk_oasis = .TRUE.  !: OASIS used
@@ -80 +81 @@
    INTEGER , PUBLIC, PARAMETER ::   jp_mfs     = 6        !: MFS  bulk                     formulation
    INTEGER , PUBLIC, PARAMETER ::   jp_none    = 7        !: for OPA when doing coupling via SAS module
+   INTEGER , PUBLIC, PARAMETER ::   jp_sbcsurge= 8        !: for OPA when doing coupling via SAS module
    INTEGER , PUBLIC, PARAMETER ::   jp_esopa   = -1       !: esopa test, ALL formulations
    

branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

@@ -91 +91 @@
    REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements
    REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements
-   LOGICAL  ::   ln_charnock ! logical flag for charnock wind stress in surge model(true) or not(false)
 
    !! * Substitutions
@@ -152 +151 @@
          &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
          &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           &
-         &                  sn_tdif, rn_zqt , rn_zu, ln_charnock 
+         &                  sn_tdif, rn_zqt , rn_zu
       !!---------------------------------------------------------------------
       !
@@ -248 +247 @@
       INTEGER  ::   ji, jj               ! dummy loop indices
       REAL(wp) ::   zcoef_qsatw, zztmp   ! local variable
-      REAL(wp) ::   z_z0, z_Cd1          ! local variable
-      REAL(wp) ::   i                    ! local variable
-      REAL(wp) ::   charn_const=0.0275    ! local variable
       REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at T-point
       REAL(wp), DIMENSION(:,:), POINTER ::   zqsatw            ! specific humidity at pst
@@ -307 +303 @@
       ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
       zztmp = 1. - albo
-#if defined key_surge
-      qsr(:,:)=0._wp
-      zqlw(:,:) = 0._wp
-#else
       IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
       ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
@@ -316 +308 @@
 
       zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
-#endif
       ! ----------------------------------------------------------------------------- !
       !     II    Turbulent FLUXES                                                    !
       ! ----------------------------------------------------------------------------- !
-      IF (ln_charnock) THEN
-          Cd(:,:)=0.0001_wp
-          DO jj = 1,jpj
-             DO ji = 1,jpi
-                z_Cd1=0._wp
-                i=1
-                !Iterate
-                DO WHILE((abs(Cd(ji,jj)-z_Cd1))>1E-6)
-                z_Cd1=Cd(ji,jj)
-                z_z0=charn_const*z_Cd1*wndm(ji,jj)**2/grav
-                Cd(ji,jj)=(0.41_wp/log(10._wp/z_z0))**2
-                i=i+1
-                ENDDO
-             ENDDO
-          ENDDO
-      ELSE
-
-        ! ... specific humidity at SST and IST
-        zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) ) 
-
-        ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point :
-        CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   &
-           &               Cd, Ch, Ce, zt_zu, zq_zu )
-
-      ENDIF
+
+      ! ... specific humidity at SST and IST
+      zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) ) 
+
+      ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point :
+      CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   &
+         &               Cd, Ch, Ce, zt_zu, zq_zu )
 
       ! ... tau module, i and j component
@@ -379 +352 @@
       !  Turbulent fluxes over ocean
       ! -----------------------------
-#if ! defined key_surge
       IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
          !! q_air and t_air are (or "are almost") given at 10m (wind reference height)
@@ -391 +363 @@
       ENDIF
       zqla (:,:) = Lv * zevap(:,:)                                                              ! Latent Heat
-#endif
 
       IF(ln_ctl) THEN
@@ -408 +379 @@
       ! ----------------------------------------------------------------------------- !
       !
-#if defined key_surge
-      emp (:,:) = 0._wp
-      qns(:,:)  = 0._wp
-#else
       emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
          &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
@@ -422 +389 @@
          &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
          &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1)
-#endif
       !
 #if defined key_lim3

branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

@@ -34 +34 @@
    USE sbcblk_core      ! surface boundary condition: bulk formulation : CORE
    USE sbcblk_mfs       ! surface boundary condition: bulk formulation : MFS
+   USE sbcsurge         ! surface boundary condition: charnock forumaltion for surge model
    USE sbcice_if        ! surface boundary condition: ice-if sea-ice model
    USE sbcice_lim       ! surface boundary condition: LIM 3.0 sea-ice model
@@ -86 +87 @@
       INTEGER ::   icpt   ! local integer
       !!
-      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   &
+      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_sbcsurge, ln_mixcpl,   &
          &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   &
          &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   &
@@ -130 +131 @@
          WRITE(numout,*) '              CORE bulk  formulation                     ln_blk_core = ', ln_blk_core
          WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs
+         WRITE(numout,*) '              surge      formulation                     ln_sbcsurge = ', ln_sbcsurge
          WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl      = ', ln_cpl
          WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl   = ', ln_mixcpl
@@ -174 +176 @@
       IF( ln_mixcpl .AND. nn_components /= jp_iam_nemo )    &
          &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) is not yet working with sas-opa coupling via oasis' )
+      IF ( lk_surge .and. .not. ( ln_sbcsurge .or. ln_ana ) ) &
+         &   CALL ctl_stop( ' surge model only compatible with analytical fluxes or surge formulae' )     
 
       !                              ! allocate sbc arrays
@@ -237 +241 @@
       IF( ln_blk_core     ) THEN   ;   nsbc = jp_core    ; icpt = icpt + 1   ;   ENDIF       ! CORE bulk            formulation
       IF( ln_blk_mfs      ) THEN   ;   nsbc = jp_mfs     ; icpt = icpt + 1   ;   ENDIF       ! MFS  bulk            formulation
+      IF( ln_sbcsurge     ) THEN   ;   nsbc = jp_sbcsurge; icpt = icpt + 1   ;   ENDIF       ! MFS  bulk            formulation
       IF( ll_purecpl      ) THEN   ;   nsbc = jp_purecpl ; icpt = icpt + 1   ;   ENDIF       ! Pure Coupled         formulation
       IF( cp_cfg == 'gyre') THEN   ;   nsbc = jp_gyre                        ;   ENDIF       ! GYRE analytical      formulation
@@ -250 +255 @@
          nstop = nstop + 1
       ENDIF
-      IF ( lk_surge .and. .not. ( ln_blk_core .or. ln_ana ) ) &
-         &   CALL ctl_stop( ' surge model only compatible with analytical fluxes or core formulae' )     
       IF(lwp) THEN
          WRITE(numout,*)
@@ -262 +265 @@
          IF( nsbc == jp_purecpl )   WRITE(numout,*) '              pure coupled formulation'
          IF( nsbc == jp_mfs     )   WRITE(numout,*) '              MFS Bulk formulation'
+         IF( nsbc == jp_sbcsurge)   WRITE(numout,*) '              Surge formulation'
          IF( nsbc == jp_none    )   WRITE(numout,*) '              OPA coupled to SAS via oasis'
          IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed formulation'
@@ -372 +376 @@
                                                                         !
       CASE( jp_mfs   )   ;   CALL sbc_blk_mfs ( kt )                    ! bulk formulation : MFS for the ocean
+      CASE( jp_sbcsurge );   CALL sbc_surge   ( kt )                    ! surge formulation
       CASE( jp_none  ) 
          IF( nn_components == jp_iam_opa ) &