Changeset 109

Timestamp:

08/07/12 16:57:48 (12 years ago)

Author:

dubos

Message:

Added Rayleigh friction
Tested : DCMIP 2.1 and DCMIP2.2 (with/without)

Location:

codes/icosagcm/trunk/src

Files:

• dissip_gcm.f90 (modified) (8 diffs)
• timeloop_gcm.f90 (modified) (1 diff)

codes/icosagcm/trunk/src/dissip_gcm.f90

@@ -7 +7 @@
   TYPE(t_field),POINTER,SAVE :: f_due_diss2(:)
 
-  TYPE(t_field),POINTER,SAVE :: f_theta(:)
+  TYPE(t_field),POINTER,SAVE :: f_theta(:), f_phi(:), f_pk(:), f_pks(:), f_p(:)
   TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:)
   TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:)
@@ -23 +23 @@
   REAL,SAVE :: cgradrot
   REAL,SAVE :: cdivgrad
+
+  INTEGER, SAVE :: rayleigh_friction_type, rayleigh_shear
+  REAL, save    :: rayleigh_tau
   
   INTEGER,SAVE :: idissip
@@ -38 +41 @@
     CALL allocate_field(f_dtheta_diss,field_t,type_real,llm)
     CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm)
+
+    CALL allocate_field(f_phi,field_t,type_real,llm)
+    CALL allocate_field(f_pk,field_t,type_real,llm)
+    CALL allocate_field(f_p,field_t,type_real,llm+1)
+    CALL allocate_field(f_pks,field_t,type_real)
     
     ALLOCATE(tau_graddiv(llm))
@@ -66 +74 @@
   REAL(rstd)            :: zz, zvert, fact
   INTEGER               :: l
+  CHARACTER(len=255)    :: rayleigh_friction_key
   
             
   INTEGER :: i,j,n,ind,it,iter
+
+   rayleigh_friction_key='none'
+   CALL getin("rayleigh_friction_type",rayleigh_friction_key)
+   SELECT CASE(TRIM(rayleigh_friction_key))
+   CASE('none')
+      rayleigh_friction_type=0
+      PRINT *, 'No Rayleigh friction'
+   CASE('dcmip2_schaer_noshear')
+      rayleigh_friction_type=1
+      rayleigh_shear=0
+      PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
+   CASE('dcmip2_schaer_shear')
+      rayleigh_shear=1
+      rayleigh_friction_type=2
+      PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
+   CASE DEFAULT
+      PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip'
+      STOP
+   END SELECT
+
+   IF(rayleigh_friction_type>0) THEN
+      rayleigh_tau=0.
+      CALL getin("rayleigh_friction_tau",rayleigh_tau)
+      rayleigh_tau = rayleigh_tau / scale_factor
+      IF(rayleigh_tau<=0) THEN
+         PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau
+         STOP
+      END IF
+   END IF
 
     CALL allocate_dissip
@@ -382 +420 @@
  
   
-  SUBROUTINE dissip(f_ue,f_due,f_ps,f_theta_rhodz,f_dtheta_rhodz)
+  SUBROUTINE dissip(f_ue,f_due,f_ps,f_phis,f_theta_rhodz,f_dtheta_rhodz)
   USE icosa
   USE theta2theta_rhodz_mod
+  USE pression_mod
+  USE exner_mod
+  USE geopotential_mod
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ue(:)
     TYPE(t_field),POINTER :: f_due(:)
-    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_ps(:), f_phis(:)
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
 
     REAL(rstd),POINTER         :: due(:,:)
+    REAL(rstd),POINTER         :: phi(:,:), ue(:,:)
     REAL(rstd),POINTER         :: due_diss1(:,:)
     REAL(rstd),POINTER         :: due_diss2(:,:)
@@ -398 +440 @@
     REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:)
 
-    INTEGER :: ind
+    INTEGER :: ind, shear
     INTEGER :: l,i,j,n
     
@@ -407 +449 @@
     CALL theta2theta_rhodz(f_ps,f_dtheta_diss,f_dtheta_rhodz_diss)
     
+    IF(rayleigh_friction_type>0) THEN
+       CALL pression(f_ps, f_p)
+       CALL exner(f_ps, f_p, f_pks, f_pk)
+       CALL geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)
+    END IF
+
     DO ind=1,ndomain
       CALL swap_dimensions(ind)
@@ -430 +478 @@
         ENDDO
       ENDDO
-    ENDDO
+
+      IF(rayleigh_friction_type>0) THEN
+         phi=f_phi(ind)
+         ue=f_ue(ind)
+         DO l=1,llm
+            DO j=jj_begin,jj_end
+               DO i=ii_begin,ii_end
+                  n=(j-1)*iim+i
+                  CALL relax(t_right, u_right)
+                  CALL relax(t_lup,   u_lup)
+                  CALL relax(t_ldown, u_ldown)
+               ENDDO
+            ENDDO
+         END DO
+      END IF
+   END DO
+
+    CONTAINS
+      SUBROUTINE relax(shift_t, shift_u)
+        USE dcmip_initial_conditions_test_1_2_3
+        REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain
+             p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain
+             fz, u3d(3), uref
+        REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values
+        LOGICAL :: hybrid_eta
+        INTEGER :: shift_u, shift_t, zcoords, nn
+        z = (phi(n,l)+phi(n+shift_t,l))/(2.*g)
+        IF(z>zh) THEN  ! relax only in the sponge layer z>zh
+!           PRINT *, 'z>zh : z,zh,l',z,zh,l
+!           STOP
+           nn = n+shift_u
+           CALL xyz2lonlat(xyz_e(nn,:),lon,lat)
+           zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
+           CALL test2_schaer_mountain(lon,lat,p,z,zcoords,hybrid_eta, &
+                hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q)
+!           u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:)
+           u3d = ulon*elon_e(nn,:) ! ulat=0
+           uref = sum(u3d*ep_e(nn,:))
+
+           fz = sin((pi/2)*(z-zh)/(ztop-zh))
+           fz = fz*fz/rayleigh_tau
+!           fz = 1/rayleigh_tau
+           due(nn,l) = due(nn,l) - fz*(ue(nn,l)-uref)
+!           due(nn,l) = due(nn,l) - fz*ue(nn,l)
+         END IF
+       END SUBROUTINE relax
       
   END SUBROUTINE dissip

codes/icosagcm/trunk/src/timeloop_gcm.f90

@@ -129 +129 @@
 
       CASE ('adam_bashforth')
-        CALL dissip(f_u,f_du,f_ps,f_theta_rhodz,f_dtheta_rhodz)
+        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
         CALL adam_bashforth_scheme