Changeset 133

Timestamp:

02/11/13 12:47:36 (11 years ago)

Author:

dubos

Message:

Working on dyn/adv coupling

Location:

codes/icosagcm/trunk/src

Files:

• advect_tracer.f90 (modified) (1 diff)
• caldyn_gcm.f90 (modified) (3 diffs)
• exner.f90 (modified) (2 diffs)
• timeloop_gcm.f90 (modified) (13 diffs)

codes/icosagcm/trunk/src/advect_tracer.f90

@@ -32 +32 @@
 
   END SUBROUTINE init_advect_tracer
-
-  SUBROUTINE advect_tracer_rhodz(f_ps, f_rhodz)
-    USE icosa
-    USE advect_mod
-    USE disvert_mod
-    USE mpipara
-    IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ps(:)     ! surface pressure, IN
-    TYPE(t_field),POINTER :: f_rhodz(:)  ! mass, OUT
-    REAL(rstd),POINTER :: rhodz(:,:) 
-    REAL(rstd),POINTER :: ps(:)
-    INTEGER :: ind
-
-    DO ind=1,ndomain
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       rhodz=f_rhodz(ind)
-       ps=f_ps(ind)
-       CALL compute_rhodz(ps,rhodz)
-    END DO
-  END SUBROUTINE advect_tracer_rhodz
-
-  SUBROUTINE compute_rhodz(ps, rhodz)
-    USE icosa
-    USE disvert_mod
-    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
-    REAL(rstd), INTENT(OUT) :: rhodz(iim*jjm,llm)
-    INTEGER :: l,i,j,ij
-    DO l = 1, llm
-       DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-             ij=(j-1)*iim+i
-             rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
-          ENDDO
-       ENDDO
-    ENDDO
-  END SUBROUTINE compute_rhodz
 
   SUBROUTINE advect_tracer(f_ps,f_u,f_q)

codes/icosagcm/trunk/src/caldyn_gcm.f90

@@ -41 +41 @@
     SELECT CASE(TRIM(def))
     CASE('lmdz')
-       caldyn_exner=1
+       caldyn_exner=lmdz
     CASE('direct')
-       caldyn_exner=2
+       caldyn_exner=direct
     CASE DEFAULT
        IF (is_mpi_root) PRINT*,'Bad selector for variable caldyn_exner : <', TRIM(def),'> options are <lmdz>, <direct>'
@@ -53 +53 @@
     SELECT CASE(TRIM(def))
     CASE('lmdz')
-       caldyn_hydrostat=1
+       caldyn_hydrostat=lmdz
     CASE('direct')
-       caldyn_hydrostat=2
+       caldyn_hydrostat=direct
     CASE DEFAULT
        IF (is_mpi_root) PRINT*,'Bad selector for variable caldyn_hydrostat : <', TRIM(def),'> options are <lmdz>, <direct>'
@@ -304 +304 @@
     REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:), Ftheta(:,:)  ! mass flux, theta flux
     REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)  ! mass flux convergence
-    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)      ! vertical velocity      
+    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)      ! vertical mass flux
     REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function
     

codes/icosagcm/trunk/src/exner.f90

@@ -2 +2 @@
 
   INTEGER :: caldyn_exner
+  INTEGER, PARAMETER :: lmdz=3, direct=4
 
 CONTAINS
@@ -45 +46 @@
     REAL(rstd) :: delta 
     
-    IF(caldyn_exner == 1) THEN          ! LMD-Z style calculation of Exner pressure
+    IF(caldyn_exner == lmdz) THEN          ! LMD-Z style calculation of Exner pressure
        !! Compute Alpha and Beta
        

codes/icosagcm/trunk/src/timeloop_gcm.f90

@@ -1 +1 @@
 MODULE timeloop_gcm_mod
-  
+
+  PRIVATE
+  
+  PUBLIC :: timeloop
+
+  INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3
+
 CONTAINS
   
@@ -24 +30 @@
   TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
   TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)
+  TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)
 
   REAL(rstd),POINTER :: phis(:)
@@ -29 +36 @@
   REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
   REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
-  REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
+  REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
   REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
   REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
   REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
+  REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
 
 !  REAL(rstd) :: dt, run_length
   INTEGER :: ind
-  INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period
-  CHARACTER(len=255) :: scheme
+  INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
+  CHARACTER(len=255) :: scheme_name
 !  INTEGER :: itaumax
 !  REAL(rstd) ::write_period
@@ -70 +78 @@
   CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
   CALL allocate_field(f_rhodz,field_t,type_real,llm)
-
-  scheme='runge_kutta'
-  CALL getin('scheme',scheme)
-
-  SELECT CASE (TRIM(scheme))
+  CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
+  CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn
+  CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
+  CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
+
+  scheme_name='runge_kutta'
+  CALL getin('scheme',scheme_name)
+
+  SELECT CASE (TRIM(scheme_name))
   CASE('euler')
+     scheme=euler
      nb_stage=1
   CASE ('runge_kutta')
+     scheme=rk4
      nb_stage=4
   CASE ('leapfrog_matsuno')
+     scheme=mlf
      matsuno_period=5
      CALL getin('matsuno_period',matsuno_period)
@@ -88 +103 @@
      CALL allocate_field(f_um2,field_u,type_real,llm)
   CASE default
-     PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             &
+     PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             &
           ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
      STOP
@@ -120 +135 @@
   CALL transfert_request(f_q,req_i1) 
 
-  CALL advect_tracer_rhodz(f_ps, f_rhodz)  ! save rhodz for transport scheme before dynamics update ps
+  DO ind=1,ndomain
+     CALL swap_dimensions(ind)
+     CALL swap_geometry(ind)
+     hfluxt=f_hfluxt(ind); hfluxt=0. ! clear accumulated fluxes
+     wfluxt=f_hfluxt(ind); wfluxt=0.
+     rhodz=f_rhodz(ind); ps=f_ps(ind)
+     CALL compute_rhodz(ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
+  END DO
 
   DO it=0,itaumax
@@ -135 +157 @@
 !       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
 !            f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
-!            f_Fe, f_W, f_dps, f_dtheta_rhodz, f_du)
+!            f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
        CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
             f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
             f_dps, f_dtheta_rhodz, f_du)
-       SELECT CASE (TRIM(scheme))
-       CASE('euler')
-          CALL  euler_scheme(.TRUE.)
-          
-       CASE ('runge_kutta')
+       SELECT CASE (scheme)
+       CASE(euler)
+          CALL euler_scheme(.TRUE.)
+       CASE (rk4)
           CALL rk_scheme(stage)
-          
-       CASE ('leapfrog_matsuno')
+       CASE (mlf)
           CALL  leapfrog_matsuno_scheme(stage)
           
@@ -155 +175 @@
           !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
           !        CALL adam_bashforth_scheme
-       CASE default
-          !        PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             &
-          !               ' > options are <euler>, <runge_kutta>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>'
-          PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             &
-               ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
+       CASE DEFAULT
           STOP
-          
        END SELECT
     END DO
@@ -168 +183 @@
     CALL euler_scheme(.FALSE.)
 
-!    CALL advect_tracer(f_Fe,f_W,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
+    IF(MOD(it+1,itau_adv)==0) THEN
+!       CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
+       DO ind=1,ndomain
+          hfluxt=f_hfluxt(ind)
+          wfluxt=f_hfluxt(ind)
+          hfluxt=0.               ! clear accumulated fluxes
+          wfluxt=0.
+       END DO
+    END IF
+
 !    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
     ENDDO
@@ -183 +207 @@
           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)
        END IF
        u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
@@ -189 +216 @@
        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
     ENDDO
-      
+
     END SUBROUTINE Euler_scheme
 
@@ -212 +239 @@
          u(:,:)=um1(:,:)+tau*du(:,:)
          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
+         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)
+         END IF
       END DO
-
     END SUBROUTINE RK_scheme
 
@@ -318 +349 @@
 
   END SUBROUTINE timeloop    
+
+  SUBROUTINE compute_rhodz(ps, rhodz)
+    USE icosa
+    USE disvert_mod
+    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
+    REAL(rstd), INTENT(OUT) :: rhodz(iim*jjm,llm)
+    INTEGER :: l,i,j,ij
+    DO l = 1, llm
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+             ij=(j-1)*iim+i
+             rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
+          ENDDO
+       ENDDO
+    ENDDO
+  END SUBROUTINE compute_rhodz
+
+  SUBROUTINE accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,tau)
+    USE icosa
+    REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
+    REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
+    REAL(rstd), INTENT(IN) :: tau
+    hfluxt = hfluxt + tau*hflux
+    wfluxt = wfluxt + tau*wflux
+  END SUBROUTINE accumulate_fluxes
   
 END MODULE timeloop_gcm_mod