Changeset 252

Timestamp:

07/23/14 14:06:47 (10 years ago)

Author:

dubos

Message:

Fix recently introduced bug : slope limiter

Location:

codes/icosagcm/trunk/src

Files:

• advect.f90 (modified) (6 diffs)
• advect_tracer.f90 (modified) (6 diffs)

codes/icosagcm/trunk/src/advect.f90

@@ -11 +11 @@
   !==========================================================================
 
-  SUBROUTINE init_advect(normal,tangent,one_over_sqrt_leng)
+  SUBROUTINE init_advect(normal,tangent,sqrt_leng)
     IMPLICIT NONE
     REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3)
     REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3)
-    REAL(rstd),INTENT(OUT) :: one_over_sqrt_leng(iim*jjm)
+    REAL(rstd),INTENT(OUT) :: sqrt_leng(iim*jjm)
 
     INTEGER :: ij
@@ -40 +40 @@
           tangent(ij+u_ldown,:)=tangent(ij+u_ldown,:)/sqrt(sum(tangent(ij+u_ldown,:)**2)+1e-50)
 
-          one_over_sqrt_leng(ij) = 1./sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))**2), &
-                                          sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))**2),  &
-                                          sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))**2)) )
+          sqrt_leng(ij) = sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))**2), &
+                                   sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))**2),  &
+                                   sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))**2)) )
     ENDDO
 
@@ -49 +49 @@
   !=======================================================================================
 
-  SUBROUTINE compute_gradq3d(qi_in,one_over_sqrt_leng_in,gradq3d_out,xyz_i,xyz_v)
+  SUBROUTINE compute_gradq3d(qi_in,sqrt_leng_in,gradq3d_out,xyz_i,xyz_v)
     USE trace
     USE omp_para
     IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: qi_in(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: one_over_sqrt_leng_in(iim*jjm)
+    REAL(rstd),INTENT(IN)  :: sqrt_leng_in(iim*jjm)
     REAL(rstd),INTENT(IN)  :: xyz_i(iim*jjm,3)
     REAL(rstd),INTENT(IN)  :: xyz_v(2*iim*jjm,3)
@@ -66 +66 @@
     INTEGER :: ij,k,ind,l
     REAL(rstd)  :: qi(iim*jjm,llm)
-    REAL(rstd)  :: one_over_sqrt_leng(iim*jjm)
+    REAL(rstd)  :: sqrt_leng(iim*jjm)
     REAL(rstd) :: gradq3d(iim*jjm,llm,3) 
     REAL(rstd) :: detx,dety,detz,det
@@ -74 +74 @@
 
     qi=qi_in
-    one_over_sqrt_leng=one_over_sqrt_leng_in
+    sqrt_leng=sqrt_leng_in
     
     CALL trace_start("compute_gradq3d1")
@@ -252 +252 @@
              maggrd = gradq3d(ij,l,1)*gradq3d(ij,l,1) + gradq3d(ij,l,2)*gradq3d(ij,l,2) + gradq3d(ij,l,3)*gradq3d(ij,l,3) 
              maggrd = sqrt(maggrd) 
-             maxq_c = qi(ij,l) + maggrd*one_over_sqrt_leng(ij)
-             minq_c = qi(ij,l) - maggrd*one_over_sqrt_leng(ij)
+             maxq_c = qi(ij,l) + maggrd*sqrt_leng(ij)
+             minq_c = qi(ij,l) - maggrd*sqrt_leng(ij)
              maxq = max(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
                   qi(ij+t_rdown,l),qi(ij+t_ldown,l))

codes/icosagcm/trunk/src/advect_tracer.f90

@@ -8 +8 @@
   TYPE(t_field),SAVE,POINTER :: f_gradq3d(:)
   TYPE(t_field),SAVE,POINTER :: f_cc(:)  ! starting point of backward-trajectory (Miura approach)
-  TYPE(t_field),SAVE,POINTER :: f_one_over_sqrt_leng(:)
+  TYPE(t_field),SAVE,POINTER :: f_sqrt_leng(:)
 
   TYPE(t_message),SAVE :: req_u, req_cc, req_wfluxt, req_q, req_rhodz, req_gradq3d
@@ -28 +28 @@
     REAL(rstd),POINTER :: tangent(:,:)
     REAL(rstd),POINTER :: normal(:,:)
-    REAL(rstd),POINTER :: one_over_sqrt_leng(:)
+    REAL(rstd),POINTER :: sqrt_leng(:)
     INTEGER :: ind
 
@@ -35 +35 @@
     CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d')
     CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc')
-    CALL allocate_field(f_one_over_sqrt_leng,field_t,type_real, name='one_over_sqrt_leng')
+    CALL allocate_field(f_sqrt_leng,field_t,type_real, name='sqrt_leng')
     CALL allocate_field(f_dzqw, field_t, type_real, llm, name='dzqw')
     CALL allocate_field(f_adzqw, field_t, type_real, llm, name='adzqw')
@@ -47 +47 @@
        normal=f_normal(ind)
        tangent=f_tangent(ind)
-       one_over_sqrt_leng=f_one_over_sqrt_leng(ind)
-       CALL init_advect(normal,tangent,one_over_sqrt_leng)
+       sqrt_leng=f_sqrt_leng(ind)
+       CALL init_advect(normal,tangent,sqrt_leng)
     END DO
 
@@ -66 +66 @@
     TYPE(t_field),POINTER :: f_rhodz(:)    ! mass field at beginning of macro time step
 
-    REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), one_over_sqrt_leng(:), gradq3d(:,:,:), cc(:,:,:)
+    REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), sqrt_leng(:), gradq3d(:,:,:), cc(:,:,:)
     REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:)
     REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 
@@ -145 +145 @@
           q       = f_q(ind)
           gradq3d = f_gradq3d(ind)
-          one_over_sqrt_leng=f_one_over_sqrt_leng(ind)
-          CALL compute_gradq3d(q(:,:,k),one_over_sqrt_leng,gradq3d,xyz_i,xyz_v)
+          sqrt_leng=f_sqrt_leng(ind)
+          CALL compute_gradq3d(q(:,:,k),sqrt_leng,gradq3d,xyz_i,xyz_v)
        END DO