Changeset 3622 for LMDZ6/trunk/libf

Timestamp:

Jan 15, 2020, 11:50:43 AM (5 years ago)

Author:

oboucher

Message:

Adding a diagnostics for the 100-m wind for estimating renewables

Location:

LMDZ6/trunk/libf/phylmd

Files:

• phys_output_ctrlout_mod.F90 (modified) (1 diff)
• phys_output_write_mod.F90 (modified) (6 diffs)

LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90

@@ -277 +277 @@
   TYPE(ctrl_out), SAVE :: o_wind10m = ctrl_out((/ 1, 1, 1, 10, 10, 10, 11, 11, 11, 11/), &
     'wind10m', '10-m wind speed', 'm/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_wind100m = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'wind100m', '100-m wind speed', 'm/s', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_wind10max = ctrl_out((/ 10, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'wind10max', '10m wind speed max', 'm/s', &

LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90

@@ -39 +39 @@
          o_t2m_min_mon, o_t2m_max_mon, &
          o_q2m, o_ustar, o_u10m, o_v10m, &
-         o_wind10m, o_wind10max, o_gusts, o_sicf, &
+         o_wind10m, o_wind10max, o_wind100m, o_gusts, o_sicf, &
          o_psol, o_mass, o_qsurf, o_qsol, &
          o_precip, o_rain_fall, o_rain_con, o_ndayrain, o_plul, o_pluc, o_plun, &
@@ -387 +387 @@
     USE tracinca_mod, ONLY: config_inca
 
+    USE vertical_layers_mod, ONLY: presnivs
+
     IMPLICIT NONE
 
@@ -428 +430 @@
     REAL, DIMENSION(klev,2) :: Ahyb_mid_bounds, Bhyb_mid_bounds
     INTEGER :: ilev
+    INTEGER, SAVE :: kmax_100m
+!$OMP THREADPRIVATE(kmax_100m)
 #ifndef CPP_XIOS
     REAL :: missing_val
@@ -443 +447 @@
     LOGICAL, PARAMETER :: debug_strataer=.FALSE.
 #endif
+    REAL,DIMENSION(klon,klev) :: z, dz
+    REAL,DIMENSION(klon)      :: zrho, zt
 
     ! On calcul le nouveau tau:
@@ -458 +464 @@
     CALL wxios_set_context
 #endif
+
+    IF (.NOT.vars_defined) THEN
+      kmax_100m=1
+      DO k=1, klev-1
+        IF (presnivs(k).GT.0.97*101325.) kmax_100m = k !--finding out max level for 100 m with a good margin
+      ENDDO
+    ENDIF
 
     Ahyb_bounds(1,1) = 0.
@@ -689 +702 @@
 
        CALL histwrite_phy(o_gusts, gustiness)
+
+       IF (vars_defined) THEN
+          DO k = 1, kmax_100m                                      !--we could stop much lower
+            zrho(:) = pplay(:,k)/t_seri(:,k)/RD                    ! air density in kg/m3
+            dz(:,k) = (paprs(:,k)-paprs(:,k+1))/zrho(:)/RG         ! layer thickness in m
+            IF (k==1) THEN
+              z(:,1) = (paprs(:,1)-pplay(:,1))/zrho(:)/RG          ! altitude middle of first layer in m
+              zt(:)  = dz(:,1)                                     ! altitude top of first layer in m
+            ELSE
+              z(:,k) = zt(:) + (paprs(:,k)-pplay(:,k))/zrho(:)/RG  ! altitude middle of layer k in m
+              zt(:)  = zt(:) + dz(:,k)                             ! altitude top of layer k in m
+            ENDIF
+          ENDDO
+          zx_tmp_fi2d(:)=missing_val
+          DO k=1, kmax_100m-1                                      !--we could stop much lower
+            DO i=1,klon
+              IF (z(i,k).LT.100..AND.z(i,k+1).GE.100.) THEN
+                zx_tmp_fi2d(i)=SQRT( (u_seri(i,k)+(100.-z(i,k))/(z(i,k+1)-z(i,k))*(u_seri(i,k+1)-u_seri(i,k)))**2.0 + &
+                                     (v_seri(i,k)+(100.-z(i,k))/(z(i,k+1)-z(i,k))*(v_seri(i,k+1)-v_seri(i,k)))**2.0 )
+              ENDIF
+            ENDDO
+          ENDDO
+       ENDIF
+       CALL histwrite_phy(o_wind100m, zx_tmp_fi2d)
 
        IF (vars_defined) THEN