Changeset 14592 for NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/ABL/ablmod.F90

Timestamp:

2021-03-05T17:03:57+01:00 (3 years ago)

Author:

gsamson

Message:

sbc_phy:

• 2 new functions:
  • pres_temp: compute air pressure from air temperature and surface pressure
  • theta_exner: compute air/surface potential temperature following a dry adiabatic transformation
• 3 modified functions:
  • Ri_bulk: use potential SST instead of SST + some simplifications TBD
  • bulk_formula: pass "rhoa" variable as input argument instead of calculing it with the old "gamma_moist" function inside the routine
  • update_qnsol_tau: pass "rhoa" variable as input argument to pass it to "bulk_formula" function
• global small cleaning/cosmetics

sbcblk_algo_{ecmwf,coare3p?}: "rhoa" variable updated for CS/WL options and passed as input argument to "update_qnsol_tau" function

sbcblk & sbcabl:
add calls to "pres_temp" & "theta_exner" functions to compute air & surface potential temperatures (over both ocean and sea-ice) before computing sensible heat flux
special attention needed to check if CS/WL increment is properly applied

ticket #2632

File:

• NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/ABL/ablmod.F90 (modified) (7 diffs)

NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/ABL/ablmod.F90

@@ -104 +104 @@
       !
       REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   zwnd_i, zwnd_j
+      REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   zsspt
+      REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   ztabs, zpre
       REAL(wp), DIMENSION(1:jpi      ,2:jpka)    ::   zCF
       !
@@ -142 +144 @@
       zrough(:,:) = z0_from_Cd( ght_abl(2), pCd_du(:,:) / MAX( pwndm(:,:), 0.5_wp ) ) ! #LB: z0_from_Cd is define in sbc_phy.F90...
 
+      ! sea surface potential temperature [K]
+      zsspt(:,:) = theta_exner( psst(:,:)+rt0, pslp_dta(:,:) )
+
       !
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -187 +192 @@
                DO ji = 1,jpi  ! surface boundary condition for temperature
                   zztmp1 = psen(ji, jj)
-                  zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )
+                  zztmp2 = psen(ji, jj) * zsspt(ji, jj)
 #if defined key_si3
                   zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj)
@@ -495 +500 @@
                zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2   &
                   &  + jp_alp1_dyn * zsig    + jp_alp0_dyn
-               zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points
-                                                             ! rn_Dt = rDt_abl / nn_fsbc
+               zcff  = (1._wp-zmsk) + zmsk * zcff2 * rDt_abl   ! zcff = 1 for masked points
                zcff  = zcff * rest_eq(ji,jj)
                u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  u_abl( ji, jj, jk, nt_a )   &
@@ -518 +522 @@
             zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2   &
                &  + jp_alp1_tra * zsig    + jp_alp0_tra
-            zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points
-                                                          ! rn_Dt = rDt_abl / nn_fsbc
+            zcff  = (1._wp-zmsk) + zmsk * zcff2 * rDt_abl   ! zcff = 1 for masked points
             tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta )   &
                &                                       + zcff   * pt_dta( ji, jj, jk )
@@ -586 +589 @@
       !
       DO_2D( 1, 1, 1, 1 )
-         ztemp          =  tq_abl( ji, jj, 2, nt_a, jp_ta )
-         zhumi          =  tq_abl( ji, jj, 2, nt_a, jp_qa )
-         zcff           = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) )
-         psen( ji, jj ) =    - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )   !GS: negative sign to respect aerobulk convention
-         pevp( ji, jj ) = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) )
+         ztemp          = tq_abl( ji, jj, 2, nt_a, jp_ta )
+         zhumi          = tq_abl( ji, jj, 2, nt_a, jp_qa )
+         zpre( ji, jj ) = pres_temp( zhumi, pslp_dta(ji,jj), ght_abl(2), ptpot=ztemp, pta=ztabs( ji, jj ) )
+         zcff           = rho_air( ztabs( ji, jj ), zhumi, zpre( ji, jj ) )
+         psen( ji, jj ) =    - cp_air(zhumi) * zcff * psen(ji,jj) * ( zsspt(ji,jj) - ztemp )         !GS: negative sign to respect aerobulk convention
+         pevp( ji, jj ) = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)  - zhumi ) )
          plat( ji, jj ) = - L_vap( psst(ji,jj) ) * pevp( ji, jj )
          rhoa( ji, jj ) = zcff
       END_2D
+      !!GS: debug only; to be removed
+      CALL iom_put ( "pres_zu", zpre(:,:) ) 
+      CALL iom_put ( "tabs_zu", ztabs(:,:) ) 
 
       DO_2D( 0, 1, 0, 1 )
@@ -639 +646 @@
       !    Wind stress relative to the moving ice ( U10m - U_ice )   !
       ! ------------------------------------------------------------ !
-      DO_2D( 0, 0, 0, 0 )
-         ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   &
-            &                      * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) )
-         ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   &
-            &                      * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) )
-      END_2D
-      CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp )
-      !
-      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   &
-         &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' )
+      !DO_2D( 0, 0, 0, 0 )
+      !   ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   &
+      !      &                      * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) )
+      !   ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   &
+      !      &                      * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) )
+      !END_2D
+      !CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp )
+      !!
+      !IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   &
+      !   &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' )
+
       ! ------------------------------------------------------------ !
       !    Wind stress relative to the moving ice ( U10m - U_ice )   !