Changeset 523 for codes/icosagcm/trunk/src

Timestamp:

01/19/17 01:58:59 (7 years ago)

Author:

dubos

Message:

Fixed Rayleigh friction (DCMIP 2.1)

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(:), 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(:)
@@ -49 +48 @@
     CALL allocate_field(f_due_diss1,field_u,type_real,llm)
     CALL allocate_field(f_due_diss2,field_u,type_real,llm)
-    CALL allocate_field(f_theta,field_t,type_real,llm)
     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))
     ALLOCATE(tau_gradrot(llm))    
@@ -435 +427 @@
       mintau=MIN(mintau,tau_divgrad(l))
     ENDDO
-       
+
+    IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau)
+
     IF(mintau>0) THEN
        itau_dissip=INT(mintau/dt)
@@ -444 +438 @@
     END IF
     itau_dissip=MAX(1,itau_dissip)
-    IF (is_master) PRINT *,"mintau ",mintau,"itau_dissip",itau_dissip," dtdissip ",dtdissip
+    IF (is_master) PRINT *,"rayleigh_tau",rayleigh_tau, "mintau ",mintau, &
+         "itau_dissip",itau_dissip," dtdissip ",dtdissip
 
   END SUBROUTINE init_dissip 
  
   
-  SUBROUTINE dissip(f_ue,f_due,f_mass,f_phis,f_theta_rhodz,f_dtheta_rhodz)
+  SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due)
   USE icosa
   USE theta2theta_rhodz_mod
@@ -459 +454 @@
   USE omp_para
   IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ue(:)
-    TYPE(t_field),POINTER :: f_due(:)
-    TYPE(t_field),POINTER :: f_mass(:), f_phis(:)
-    TYPE(t_field),POINTER :: f_theta_rhodz(:)
-    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
+    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:)
+    TYPE(t_field),POINTER :: f_ue(:), f_due(:)
 
     REAL(rstd),POINTER         :: due(:,:)
@@ -483 +476 @@
     CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss)
     
-! later for openmp    
-!    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
       IF (.NOT. assigned_domain(ind)) CYCLE
@@ -515 +501 @@
 !      due=0
 
-! later for openmp  
-!      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
+      IF(rayleigh_friction_type>0) THEN
+         phi=f_geopot(ind)
+         ue=f_ue(ind)
+         DO l=1,llm
+            DO ij=ij_begin,ij_end
+               CALL relax(t_right, u_right)
+               CALL relax(t_lup,   u_lup)
+               CALL relax(t_ldown, u_ldown)
+            ENDDO
+         END DO
+      END IF
    END DO
 
@@ -558 +540 @@
            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) - itau_dissip*fz*(ue(nn,l)-uref)
+!           fz = 1./rayleigh_tau
 !           due(nn,l) = due(nn,l) - fz*ue(nn,l)
          END IF

codes/icosagcm/trunk/src/timeloop_gcm.f90

@@ -284 +284 @@
                f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
        
-          CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz)
+          CALL dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_u, f_dtheta_rhodz,f_du)
           
           CALL euler_scheme(.FALSE.)  ! update only u, theta