Changeset 138 for codes/icosagcm/trunk/src/timeloop_gcm.f90

Timestamp:

02/16/13 17:03:57 (11 years ago)

Author:

dubos

Message:

Transport now working again - tested with dcmip4.1.0

File:

• codes/icosagcm/trunk/src/timeloop_gcm.f90 (modified) (11 diffs)

codes/icosagcm/trunk/src/timeloop_gcm.f90

@@ -46 +46 @@
   INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
   CHARACTER(len=255) :: scheme_name
-  LOGICAL :: fluxt_zero ! set to .TRUE. to start accumulating fluxes in time
+  LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
 !  INTEGER :: itaumax
 !  REAL(rstd) ::write_period
@@ -141 +141 @@
      CALL swap_geometry(ind)
      rhodz=f_rhodz(ind); ps=f_ps(ind)
-     CALL compute_rhodz(ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
+     CALL compute_rhodz(.TRUE., ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
   END DO
-  fluxt_zero=.FALSE.
-
+  fluxt_zero=.TRUE.
+
+  ! check that rhodz is consistent with ps
+  CALL transfert_request(f_rhodz,req_i1)
+  CALL transfert_request(f_ps,req_i1)
+  DO ind=1,ndomain
+     CALL swap_dimensions(ind)
+     CALL swap_geometry(ind)
+     rhodz=f_rhodz(ind); ps=f_ps(ind)
+     CALL compute_rhodz(.FALSE., ps, rhodz)   
+  END DO
+  
   DO it=0,itaumax
 
@@ -182 +192 @@
 
     IF(MOD(it+1,itau_adv)==0) THEN
+       CALL transfert_request(f_wfluxt,req_i1) ! FIXME
+!       CALL transfert_request(f_hfluxt,req_e1) ! FIXME
+
        CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
        fluxt_zero=.TRUE.
+
+       ! FIXME : check that rhodz is consistent with ps
+       CALL transfert_request(f_rhodz,req_i1)
+       CALL transfert_request(f_ps,req_i1)
+       CALL transfert_request(f_dps,req_i1) ! FIXME
+       CALL transfert_request(f_wflux,req_i1) ! FIXME
+       DO ind=1,ndomain
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          rhodz=f_rhodz(ind); ps=f_ps(ind); dps=f_dps(ind); 
+          wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
+          CALL compute_rhodz(.FALSE., ps, rhodz)   
+       END DO
+
     END IF
 
@@ -195 +222 @@
     LOGICAL :: with_dps
     INTEGER :: ind
-      
     DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
        IF(with_dps) THEN
           ps=f_ps(ind) ; dps=f_dps(ind) ; 
           ps(:)=ps(:)+dt*dps(:)
           hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
-          wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind)
-          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero)
+          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
+          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
        END IF
        u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
@@ -221 +249 @@
       
       DO ind=1,ndomain
+         CALL swap_dimensions(ind)
+         CALL swap_geometry(ind)
          ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
          psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
@@ -234 +264 @@
          IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
             hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
-            wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind)
-            CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero)
+            wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
+            CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
          END IF
       END DO
@@ -245 +275 @@
 
       DO ind=1,ndomain
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
@@ -277 +309 @@
     tau = dt/nb_stage
       DO ind=1,ndomain
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
+
         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
@@ -320 +355 @@
 
       DO ind=1,ndomain
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
@@ -343 +380 @@
   END SUBROUTINE timeloop    
 
-  SUBROUTINE compute_rhodz(ps, rhodz)
+  SUBROUTINE compute_rhodz(comp, ps, rhodz)
     USE icosa
     USE disvert_mod
+    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
     REAL(rstd), INTENT(IN) :: ps(iim*jjm)
     REAL(rstd), INTENT(OUT) :: rhodz(iim*jjm,llm)
-    INTEGER :: l,i,j,ij
+    REAL(rstd) :: m, err
+    INTEGER :: l,i,j,ij,dd
+    err=0.
+    IF(comp) THEN 
+       dd=1
+    ELSE
+!       dd=-1
+       dd=0
+    END IF
+
     DO l = 1, llm
-       DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
+       DO j=jj_begin-dd,jj_end+dd
+          DO i=ii_begin-dd,ii_end+dd
              ij=(j-1)*iim+i
-             rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
+             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
+             IF(comp) THEN
+                rhodz(ij,l) = m
+             ELSE
+                err = MAX(err,abs(m-rhodz(ij,l)))
+             END IF
           ENDDO
        ENDDO
     ENDDO
+
+    IF(.NOT. comp) THEN
+       IF(err>1e-10) THEN
+          PRINT *, 'Discrepancy between ps and rhodz detected', err
+          STOP
+       ELSE
+!          PRINT *, 'No discrepancy between ps and rhodz detected'
+       END IF
+    END IF
+
   END SUBROUTINE compute_rhodz
 
-  SUBROUTINE accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,tau,fluxt_zero)
+  SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
     USE icosa
     REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
@@ -366 +428 @@
     LOGICAL, INTENT(INOUT) :: fluxt_zero
     IF(fluxt_zero) THEN
+!       PRINT *, 'Accumulating fluxes (first)'
        fluxt_zero=.FALSE.
        hfluxt = tau*hflux
        wfluxt = tau*wflux
     ELSE
+!       PRINT *, 'Accumulating fluxes (next)'
        hfluxt = hfluxt + tau*hflux
        wfluxt = wfluxt + tau*wflux