Changeset 192 for codes/icosagcm/trunk

Timestamp:

06/25/14 15:00:24 (10 years ago)

Author:

dubos

Message:

Fixed DCMIP5 physics/etat0

Location:

codes/icosagcm/trunk/src

Files:

• etat0.f90 (modified) (2 diffs)
• etat0_dcmip5.f90 (modified) (1 diff)
• physics_dcmip.f90 (modified) (12 diffs)

codes/icosagcm/trunk/src/etat0.f90

@@ -7 +7 @@
   SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
     USE icosa
+    USE mpipara, ONLY : is_mpi_root
     USE disvert_mod
     USE etat0_williamson_mod, ONLY : etat0_williamson_new
@@ -61 +62 @@
        CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
      CASE ('dcmip4')
+        IF(nqtot<2) THEN
+           IF (is_mpi_root)  THEN
+              PRINT *, "nqtot must be at least 2 for test case DCMIP4"
+           END IF
+           STOP
+        END IF
        CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
      CASE ('dcmip5')

codes/icosagcm/trunk/src/etat0_dcmip5.f90

@@ -96 +96 @@
     lonc=Pi 
 
-    Tv0=T0/(1+0.608*q0)
+    Tv0=T0*(1+0.608*q0)
     Tvt=Tv0-Gamma*zt
   

codes/icosagcm/trunk/src/physics_dcmip.f90

@@ -16 +16 @@
   IMPLICIT NONE
 
-    testcase=1
+    testcase=1 ! OK for 4.2 (moist baroclinic instability)
     CALL getin("dcmip_physics",testcase)
     CALL allocate_field(f_out_i,field_t,type_real,llm)
@@ -106 +106 @@
     INTEGER :: i,j,l,ij
      
-    PRINT *,'Entering in DCMIP physics'    
+!    PRINT *,'Entering in DCMIP physics'    
     CALL compute_pression(ps,p,0)
     CALL compute_exner(ps,p,pks,pk,0)
@@ -112 +112 @@
     CALL compute_geopotential(phis,pks,pk,theta,phi,0)
     CALL compute_theta_rhodz2temperature(ps,theta_rhodz,T,0)
-   CALL compute_wind_centered(ue,uc)
-   CALL compute_wind_centered_lonlat_compound(uc, u, v)
-
+    ! Reconstruct wind vector at hexagons
+    CALL compute_wind_centered(ue,uc)
+    CALL compute_wind_centered_lonlat_compound(uc, u, v)
+
+    ! Convert from Tv to T
+    DO l=1,llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          T(ij,l)=T(ij,l)/(1+0.608*q(ij,l))
+        ENDDO
+      ENDDO
+    ENDDO       
+    
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
@@ -122 +133 @@
     ENDDO
    
+    ! bottom-up indexing (DYNAMICO) : u,utemp, v,vtemp                                                                     
+    ! top-down vertical indexing (DCMIP) : ufi, vfi, ...                                                             
+    ! => copy fields and mirror vertical index                                                       
     DO l=1,llm+1
       DO j=jj_begin,jj_end
@@ -131 +145 @@
     ENDDO
     
+    ! Pressure inside layers                                                                                             
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -140 +155 @@
     ENDDO
     
+    ! Pressure difference between two layers
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -149 +165 @@
     ENDDO       
        
-    
-!    ufi=u
-!    vfi=v    
-    
-    DO l=1,llm
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          T(ij,l)=T(ij,l)/(1+0.608*q(ij,l))
-        ENDDO
-      ENDDO
-    ENDDO       
-    
+    ! Copy T,u,v,q for input to physics
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -174 +178 @@
     ENDDO       
     
-!    q=0
-!    out_i=T
-    
     CALL simple_physics(iim*jjm, llm, dt, lat, tfi, qfi , ufi, vfi, pmid, pint, pdel, 1/pdel, ps, precl, testcase) 
     
+    ! Copy back T,u,v,q and mirror vertical index
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -192 +194 @@
 
 
+    ! Convert back T to Tv
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -201 +204 @@
     ENDDO       
 
-!    out_i=q
-        
+    ! Compute velocity update at hexagons                                                                                
     utemp=utemp-u
     vtemp=vtemp-v
 
+    ! lon-lat -> 3D
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -215 +218 @@
     ENDDO
 
-!    out_i=ufi
-    
+    ! Update velocity at velocity points
     DO l=1,llm
       DO j=jj_begin,jj_end
@@ -519 +521 @@
       end do
 
-      IF (test==0) return
 !===============================================================================
 ! Send variables to history file - THIS PROCESS WILL BE MODEL SPECIFIC