Changeset 968 for codes

Timestamp:

09/02/19 12:22:39 (5 years ago)

Author:

jisesh

Message:

devel : nudge also surface pressure

Location:

codes/icosagcm/devel/src

Files:

• dissip/guided_mod.f90 (modified) (3 diffs)
• dissip/nudging_mod.f90 (modified) (15 diffs)
• time/timeloop_gcm.f90 (modified) (1 diff)

codes/icosagcm/devel/src/dissip/guided_mod.f90

@@ -8 +8 @@
 
 
-  SUBROUTINE init_guided(f_u,f_theta_rhodz)
+  SUBROUTINE init_guided(f_u,f_theta_rhodz,f_ps)
     USE icosa
     USE guided_ncar_mod, ONLY : init_guided_ncar => init_guided
@@ -15 +15 @@
     TYPE(t_field),POINTER :: f_u(:)
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_ps(:)
     
     guided_type='none'
@@ -26 +27 @@
         
       CASE ('nudging')
-        CALL init_guided_nudging(f_u,f_theta_rhodz)
-        
+        CALL init_guided_nudging(f_u,f_theta_rhodz,f_ps)
       CASE DEFAULT
          PRINT*,"Bad selector for varaible guided_type >",TRIM(guided_type),"> option are <none>, <dcmip1>, <nudging>"

codes/icosagcm/devel/src/dissip/nudging_mod.f90

@@ -12 +12 @@
   PRIVATE
   ! nudging will be active outside a disc of radius 'radius' centered at 'center_lon', 'center lat'.
-  REAL(rstd) :: center_lon, center_lat, nudging_radius, time
-  !$OMP THREADPRIVATE(center_lon, center_lat, nudging_radius, time)
-
+  REAL(rstd) :: center_lon, center_lat, nudging_radius, nudging_relaxation_time
+  !$OMP THREADPRIVATE(center_lon, center_lat, nudging_radius, nudging_relaxation_time)
   TYPE(t_field),POINTER :: f_relax_coef_e(:), f_target_ue(:), &
-       f_relax_coef_i(:), f_target_theta_rhodz(:)
+       f_relax_coef_i(:), f_target_theta_rhodz(:),f_target_ps(:)
   CHARACTER(LEN=255),SAVE :: guided_nudging_field
   INTEGER,SAVE :: nudging_time
@@ -24 +23 @@
 CONTAINS
 
-  SUBROUTINE init_guided(f_u,f_theta_rhodz)
+  SUBROUTINE init_guided(f_u,f_theta_rhodz,f_ps)
     USE getin_mod, ONLY : getin
     USE math_const, ONLY : pi
@@ -30 +29 @@
     TYPE(t_field),POINTER :: f_u(:)! initial condition
     TYPE(t_field),POINTER :: f_theta_rhodz(:)! initial condition
+    TYPE(t_field),POINTER :: f_ps(:)! initial condition
     REAL(rstd), POINTER :: ue(:,:), target_ue(:,:), coef_e(:) 
     REAL(rstd), POINTER :: theta_rhodz(:,:,:), target_theta_rhodz(:,:,:), coef_i(:) 
+    REAL(rstd), POINTER :: ps2(:), target_ps(:)
     INTEGER :: ind
     ! read DEF keys describing how to relax
@@ -41 +42 @@
     CALL getin('nudging_radius', nudging_radius)
     nudging_radius = nudging_radius / scale_factor
-
-    nudging_time=0.
+    ! we should check that radius>0
+
+    nudging_time=0
     CALL getin('nudging_time', nudging_time)
     nudging_time = nudging_time/scale_factor
-    ! we should check that radius>0
+    !nudging_time=0
+
+    nudging_relaxation_time = 0.
+    CALL getin('nudging_relaxation_time', nudging_relaxation_time)
+    nudging_relaxation_time = nudging_relaxation_time/scale_factor
+
     CALL getin("guided_nudging_field",guided_nudging_field)
 
     SELECT CASE(TRIM(guided_nudging_field))
-        CASE ('wind')
-            CALL allocate_field(f_relax_coef_e, field_u, type_real, name='nudging_coef_e')
-            CALL allocate_field(f_target_ue, field_u, type_real, llm, name='nudging_target_e')
-        CASE ('temperature')
-            CALL allocate_field(f_relax_coef_i, field_t, type_real, name='nudging_coef_i')
-            CALL allocate_field(f_target_theta_rhodz, field_t, type_real, llm,nqdyn, name='nudging_target_theta')
-        CASE ('wind_temperature')
+        CASE('wind','temperature','wind_temperature')
             CALL allocate_field(f_relax_coef_e, field_u, type_real, name='nudging_coef_e')
             CALL allocate_field(f_target_ue, field_u, type_real, llm, name='nudging_target_e')
             CALL allocate_field(f_relax_coef_i, field_t, type_real, name='nudging_coef_i')
             CALL allocate_field(f_target_theta_rhodz, field_t, type_real, llm,nqdyn, name='nudging_target_theta')
+            CALL allocate_field(f_target_ps, field_t, type_real, name='nudging_target_ps')
         CASE DEFAULT
                  PRINT*,"Bad selector for varaible init_guided_nudging>",TRIM(guided_nudging_field)
@@ -65 +67 @@
                  STOP
     END SELECT
-
     ! compute relax_coef and target_ue
     center_lon = center_lon * pi/180.
@@ -79 +80 @@
             target_theta_rhodz = f_target_theta_rhodz(ind)
             theta_rhodz = f_theta_rhodz(ind)
+            target_ps = f_target_ps(ind)
+            ps2 = f_ps(ind)
             CALL compute_relax_coef(coef_e, coef_i)
             CALL compute_target_u(ue, target_ue)
             CALL compute_target_center(theta_rhodz, target_theta_rhodz)
+            CALL compute_target_center2(ps2, target_ps)
     END DO
   END SUBROUTINE init_guided
@@ -101 +105 @@
   FUNCTION relax_coef(lon,lat)
     USE spherical_geom_mod, ONLY : dist_lonlat
+    USE time_mod!, ONLY : dt ! time step for a full RK step
     REAL(rstd), INTENT(IN) :: lon,lat
     REAL(rstd) :: relax_coef, dist, c
     ! NB : dist is computed on unit sphere
-    CALL dist_lonlat(lon, lat, center_lon, center_lat, dist) 
-    c = tanh((1.-radius*dist/nudging_radius)*20.) ! 1 inside circle, -1 outside
-    IF (c>0.99) c=1
-    IF (c<-0.99) c=-1
-    relax_coef = .5*(1.+c) ! rescale to [0,1] range ; ! 1 inside circle, 0 outside
+    IF(nudging_radius>0.) THEN
+       CALL dist_lonlat(lon, lat, center_lon, center_lat, dist) 
+       c = tanh((1.-radius*dist/nudging_radius)*20.) ! 1 inside circle, -1 outside
+       IF (c>0.99) c=1
+       IF (c<-0.99) c=-1
+       relax_coef = .5*(1.+c) ! rescale to [0,1] range ; ! 1 inside circle, 0 outside
+    ELSE
+       c = 0.
+    END IF
+
+    relax_coef = 1. - (1.-relax_coef) * dt / (dt + nudging_relaxation_time)
   END FUNCTION relax_coef
 
@@ -139 +150 @@
   END SUBROUTINE compute_target_center
 
+  SUBROUTINE compute_target_center2(theta_rhodz, target_theta_rhodz)
+    REAL(rstd), INTENT(OUT)  :: target_theta_rhodz(iim*jjm,nqdyn)
+    REAL(rstd), INTENT(IN) :: theta_rhodz(iim*jjm,nqdyn)
+    INTEGER :: ij, iq
+    DO iq=1, nqdyn
+          DO ij=ij_begin_ext, ij_end_ext
+             target_theta_rhodz(ij,iq)=theta_rhodz(ij,iq)
+          END DO
+    END DO
+  END SUBROUTINE compute_target_center2
+
 !----------------------------- Relax towards target ------------------------------
   
@@ -146 +168 @@
     USE wind_mod
     USE transfert_mod
-    USE time_mod
+    USE time_mod, ONLY : dt ! time step for a full RK step
     USE vertical_remap_mod
     USE compute_pression_mod, ONLY : pression_mid
@@ -157 +179 @@
     TYPE(t_field),POINTER :: f_q(:)
     TYPE(t_field),POINTER, SAVE :: f_T_guided(:),  f_ulon_guided(:), f_ulat_guided(:),f_T_guided_interp(:), &
-                                   f_ulon_guided_interp(:),f_ulat_guided_interp(:)
+                                   f_ulon_guided_interp(:),f_ulat_guided_interp(:),f_ps_guided(:)
     REAL(rstd), POINTER :: target_ue(:,:), ue(:,:), coef_e(:)
     REAL(rstd), POINTER :: target_theta_rhodz(:,:,:), theta_rhodz(:,:,:), coef_i(:) 
+    REAL(rstd), POINTER :: target_ps(:), ps2(:)
     INTEGER :: ind
     
-    IF (abs(MOD(tt,REAL(nudging_time))-dt) < 1.0D-2) THEN 
-       
+    IF (abs(MOD(tt,REAL(nudging_time))-dt) < 1.0D-2) THEN
+
+      CALL allocate_field(f_ps_guided, field_t, type_real, name='nudging_ps')
+      CALL xios_read_field("ps_guided_read",f_ps_guided)
+      !f_ps = f_ps_guided  !================================================>FIX
+      CALL transfert_request(f_ps_guided,req_i0)
+      CALL transfert_request(f_ps_guided,req_i1)
       CALL allocate_field(f_T_guided, field_t, type_real, llm, name='nudging_T')
       CALL allocate_field(f_T_guided_interp, field_t, type_real, llm, name='nudging_T')
@@ -172 +200 @@
       CALL allocate_field(f_pmid_target,field_t,type_real,llm,name='target_pressure') 
       CALL xios_read_field("T_guided_read",f_T_guided)
-      CALL pression_mid(f_ps, f_pmid_target)
+      CALL pression_mid(f_ps, f_pmid_target) ! call ps nudging before this
       CALL vertical_remap_extdata(f_T_guided,f_pmid_target,f_T_guided_interp)
       CALL transfert_request(f_T_guided,req_i0)
@@ -191 +219 @@
       CALL transfert_request(f_ulat_guided,req_i1)
       CALL transfert_request(f_ulat_guided_interp,req_i1)
-      CALL xios_write_field("ulat_guided_out",f_ulat_guided)
+      !CALL xios_write_field("ulat_guided_out",f_ulat_guided)
+      CALL xios_write_field("ulat_guided_out",f_T_guided_interp)
       CALL xios_write_field("ulon_guided_out",f_ulon_guided)
       CALL xios_write_field("T_guided_out",f_T_guided)
+      CALL xios_write_field("PS_guided_out",f_ps_guided)
       CALL ulonlat2un(f_ulon_guided_interp, f_ulat_guided,f_target_ue) 
       CALL deallocate_field(f_T_guided)
@@ -202 +232 @@
       CALL deallocate_field(f_ulat_guided_interp)
       CALL deallocate_field(f_pmid_target)
-    
+      !CALL deallocate_field(f_ps_guided)
+      print*,"-----------------------nudging end------------------"
+
     ENDIF
 
@@ -217 +249 @@
        target_theta_rhodz = f_target_theta_rhodz(ind)
        theta_rhodz = f_theta_rhodz(ind)
-       CALL compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
+       !CALL compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz)
+       target_ps = f_ps_guided(ind)
+       ps2 = f_ps(ind)
+       CALL compute_guided_center2(coef_i, target_ps, ps2)
     END DO
 
@@ -253 +288 @@
     END DO
   END SUBROUTINE compute_guided_center
+
+  SUBROUTINE compute_guided_center2(coef_i, target_theta_rhodz, theta_rhodz)
+    REAL(rstd), INTENT(IN) :: coef_i(iim*jjm)
+    REAL(rstd), INTENT(IN) :: target_theta_rhodz(iim*jjm,nqdyn)
+    REAL(rstd), INTENT(INOUT) :: theta_rhodz(iim*jjm,nqdyn)
+    INTEGER :: ij, iq
+    DO iq=1, nqdyn
+          DO ij=ij_begin_ext, ij_end_ext
+             theta_rhodz(ij,iq) = theta_rhodz(ij,iq)*coef_i(ij) + &
+                  target_theta_rhodz(ij,iq)*(1.-coef_i(ij))
+          END DO
+    END DO
+  END SUBROUTINE compute_guided_center2
   
 END MODULE nudging_mod

codes/icosagcm/devel/src/time/timeloop_gcm.f90

@@ -167 +167 @@
 
     CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
-    
-    CALL init_guided(f_u,f_theta_rhodz)
+   
+    CALL init_guided(f_u,f_theta_rhodz,f_ps)
 
     CALL transfert_request(f_phis,req_i0)