Changeset 129


Ignore:
Timestamp:
02/04/13 15:15:35 (11 years ago)
Author:
dubos
Message:

Changed timeloop design for multistage schemes

Location:
codes/icosagcm/trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/caldyn.f90

    r110 r129  
    3030  END SUBROUTINE init_caldyn 
    3131   
    32   SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
     32  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
    3333  USE icosa 
    3434  USE caldyn_gcm_mod, ONLY : caldyn_gcm=>caldyn 
    3535  USE caldyn_adv_mod, ONLY : caldyn_adv=>caldyn 
    3636  IMPLICIT NONE 
    37   INTEGER,INTENT(IN)    :: it 
     37  LOGICAL,INTENT(IN)    :: write_out 
    3838  TYPE(t_field),POINTER :: f_phis(:) 
    3939  TYPE(t_field),POINTER :: f_ps(:) 
    … …  
    4747    SELECT CASE (TRIM(caldyn_type)) 
    4848      CASE('gcm') 
    49         CALL caldyn_gcm(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
     49        CALL caldyn_gcm(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
    5050      CASE('adv') 
    51         CALL caldyn_adv(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
     51        CALL caldyn_adv(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
    5252    END SELECT 
    5353   

  • codes/icosagcm/trunk/src/caldyn_adv.f90

    r110 r129  
    7979   
    8080   
    81   SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
     81  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
    8282  USE icosa 
    8383  USE vorticity_mod 
    … …  
    8585  USE theta2theta_rhodz_mod 
    8686  IMPLICIT NONE 
    87   INTEGER,INTENT(IN)    :: it 
     87  LOGICAL,INTENT(IN)    :: write_out 
    8888  TYPE(t_field),POINTER :: f_phis(:) 
    8989  TYPE(t_field),POINTER :: f_ps(:) 
    … …  
    133133 
    134134 
    135     IF (mod(it,itau_out)==0 ) THEN 
     135    IF (write_out) THEN 
    136136      CALL writefield("ps",f_ps) 
    137137    ENDIF 

  • codes/icosagcm/trunk/src/caldyn_gcm.f90

    r128 r129  
    8080  END SUBROUTINE allocate_caldyn 
    8181    
    82   SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
     82  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
    8383    USE icosa 
    8484    USE vorticity_mod 
    … …  
    8686    USE theta2theta_rhodz_mod 
    8787    IMPLICIT NONE 
    88     INTEGER,INTENT(IN)    :: it 
     88    LOGICAL,INTENT(IN)    :: write_out 
    8989    TYPE(t_field),POINTER :: f_phis(:) 
    9090    TYPE(t_field),POINTER :: f_ps(:) 
    … …  
    168168    END SELECT 
    169169 
    170     IF (mod(it,itau_out)==0 ) THEN 
     170    IF (write_out) THEN 
    171171       PRINT *,'CALL write_output_fields' 
    172172       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, & 

  • codes/icosagcm/trunk/src/dissip_gcm.f90

    r109 r129  
    250250    ENDDO  
    251251    PRINT *,"gradiv : dumax",dumax 
    252      
     252    PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., & 
     253         'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000 
     254    PRINT *, 'Max. time step assuming c=340 m/s and Courant number=2.8 :', & 
     255         2.8/340.*dumax**(-.5/nitergdiv) 
     256 
    253257    cgraddiv=dumax**(-1./nitergdiv) 
    254258    PRINT *,"cgraddiv : ",cgraddiv 
    … …  
    469473            n=(j-1)*iim+i 
    470474 
    471             due(n+u_right,l)=due(n+u_right,l)-0.5*(due_diss1(n+u_right,l)/tau_graddiv(l)+due_diss2(n+u_right,l)/tau_gradrot(l)) 
    472             due(n+u_lup,l)=due(n+u_lup,l)    -0.5*(due_diss1(n+u_lup,l)  /tau_graddiv(l)+due_diss2(n+u_lup,l)  /tau_gradrot(l)) 
    473             due(n+u_ldown,l)=due(n+u_ldown,l)-0.5*(due_diss1(n+u_ldown,l)/tau_graddiv(l)+due_diss2(n+u_ldown,l)/tau_gradrot(l)) 
    474  
    475             dtheta_rhodz(n,l)=dtheta_rhodz(n,l)-0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l) 
    476  
     475            due(n+u_right,l) = -0.5*( due_diss1(n+u_right,l)/tau_graddiv(l) + due_diss2(n+u_right,l)/tau_gradrot(l) ) 
     476            due(n+u_lup,l)   = -0.5*( due_diss1(n+u_lup,l)  /tau_graddiv(l) + due_diss2(n+u_lup,l)  /tau_gradrot(l) ) 
     477            due(n+u_ldown,l) = -0.5*( due_diss1(n+u_ldown,l)/tau_graddiv(l) + due_diss2(n+u_ldown,l)/tau_gradrot(l) ) 
     478 
     479            dtheta_rhodz(n,l) = -0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l) 
    477480          ENDDO 
    478481        ENDDO 

  • codes/icosagcm/trunk/src/timeloop_gcm.f90

    r124 r129  
    1414  IMPLICIT NONE 
    1515  TYPE(t_field),POINTER :: f_phis(:) 
    16   TYPE(t_field),POINTER :: f_theta(:) 
     16!  TYPE(t_field),POINTER :: f_theta(:) 
    1717  TYPE(t_field),POINTER :: f_q(:) 
    1818  TYPE(t_field),POINTER :: f_dtheta(:) 
    … …  
    3535!  REAL(rstd) :: dt, run_length 
    3636  INTEGER :: ind 
    37   INTEGER :: it,i,j,n 
     37  INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period 
    3838  CHARACTER(len=255) :: scheme 
    39   INTEGER :: matsuno_period 
    4039!  INTEGER :: itaumax 
    4140!  REAL(rstd) ::write_period 
    … …  
    5453!  dt=dt/scale_factor 
    5554 
    56   scheme='adam_bashforth' 
     55! Trends 
     56  CALL allocate_field(f_dps,field_t,type_real) 
     57  CALL allocate_field(f_du,field_u,type_real,llm) 
     58  CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm) 
     59! Model state at current time step (RK/MLF/Euler) 
     60  CALL allocate_field(f_phis,field_t,type_real) 
     61  CALL allocate_field(f_ps,field_t,type_real) 
     62  CALL allocate_field(f_u,field_u,type_real,llm) 
     63  CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) 
     64! Model state at previous time step (RK/MLF) 
     65  CALL allocate_field(f_psm1,field_t,type_real) 
     66  CALL allocate_field(f_um1,field_u,type_real,llm) 
     67  CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm) 
     68! Tracers 
     69  CALL allocate_field(f_q,field_t,type_real,llm,nqtot) 
     70 
     71  scheme='runge_kutta' 
    5772  CALL getin('scheme',scheme) 
    58    
    59   matsuno_period=5 
    60   CALL getin('matsuno_period',matsuno_period) 
    61   IF (TRIM(scheme)=='leapfrog') matsuno_period=itaumax+1 
    62  
     73 
     74  SELECT CASE (TRIM(scheme)) 
     75  CASE('euler') 
     76     nb_stage=1 
     77  CASE ('runge_kutta') 
     78     nb_stage=4 
     79  CASE ('leapfrog_matsuno') 
     80     matsuno_period=5 
     81     CALL getin('matsuno_period',matsuno_period) 
     82     nb_stage=matsuno_period+1 
     83     ! Model state 2 time steps ago (MLF) 
     84     CALL allocate_field(f_psm2,field_t,type_real) 
     85     CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm) 
     86     CALL allocate_field(f_um2,field_u,type_real,llm) 
     87  CASE default 
     88     PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             & 
     89          ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>' 
     90     STOP 
     91  END SELECT 
     92   
    6393!  write_period=0 
    6494!  CALL getin('write_period',write_period) 
    … …  
    6797!  PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out 
    6898  
    69   CALL allocate_field(f_phis,field_t,type_real) 
    70    
    71   CALL allocate_field(f_ps,field_t,type_real) 
    72   CALL allocate_field(f_psm1,field_t,type_real) 
    73   CALL allocate_field(f_psm2,field_t,type_real) 
    74   CALL allocate_field(f_dps,field_t,type_real) 
    75   CALL allocate_field(f_dpsm1,field_t,type_real) 
    76   CALL allocate_field(f_dpsm2,field_t,type_real) 
    77  
    78   CALL allocate_field(f_u,field_u,type_real,llm) 
    79   CALL allocate_field(f_um1,field_u,type_real,llm) 
    80   CALL allocate_field(f_um2,field_u,type_real,llm) 
    81   CALL allocate_field(f_du,field_u,type_real,llm) 
    82   CALL allocate_field(f_dum1,field_u,type_real,llm) 
    83   CALL allocate_field(f_dum2,field_u,type_real,llm) 
    84  
    85   CALL allocate_field(f_theta,field_t,type_real,llm) 
    86   CALL allocate_field(f_dtheta,field_t,type_real,llm) 
    87  
    88   CALL allocate_field(f_q,field_t,type_real,llm,nqtot) 
    89  
    90   CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) 
    91   CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm) 
    92   CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm) 
    93   CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm) 
    94   CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm) 
    95   CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm) 
     99! Trends at previous time steps needed only by Adams-Bashforth 
     100!  CALL allocate_field(f_dpsm1,field_t,type_real) 
     101!  CALL allocate_field(f_dpsm2,field_t,type_real) 
     102!  CALL allocate_field(f_dum1,field_u,type_real,llm) 
     103!  CALL allocate_field(f_dum2,field_u,type_real,llm) 
     104!  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm) 
     105!  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm) 
     106 
     107!  CALL allocate_field(f_theta,field_t,type_real,llm) 
     108!  CALL allocate_field(f_dtheta,field_t,type_real,llm) 
    96109 
    97110  CALL init_dissip 
    … …  
    113126    ENDIF 
    114127     
    115     CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 
    116     CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
     128!    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 
     129 
     130    DO stage=1,nb_stage 
     131       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), & 
     132            f_phis,f_ps,f_theta_rhodz,f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 
     133       SELECT CASE (TRIM(scheme)) 
     134       CASE('euler') 
     135          CALL  euler_scheme 
     136           
     137       CASE ('runge_kutta') 
     138          CALL rk_scheme(stage) 
     139           
     140       CASE ('leapfrog_matsuno') 
     141          CALL  leapfrog_matsuno_scheme(stage) 
     142           
     143          !      CASE ('leapfrog') 
     144          !        CALL leapfrog_scheme 
     145          ! 
     146          !      CASE ('adam_bashforth') 
     147          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) 
     148          !        CALL adam_bashforth_scheme 
     149       CASE default 
     150          !        PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             & 
     151          !               ' > options are <euler>, <runge_kutta>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>' 
     152          PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             & 
     153               ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>' 
     154          STOP 
     155           
     156       END SELECT 
     157    END DO 
    117158!    CALL advect_tracer(f_ps,f_u,f_q) 
    118      
    119     SELECT CASE (TRIM(scheme)) 
    120       CASE('euler') 
    121         CALL  euler_scheme 
    122  
    123       CASE ('leapfrog') 
    124         CALL leapfrog_scheme 
    125  
    126       CASE ('runge_kutta') 
    127         CALL rk_scheme 
    128  
    129       CASE ('leapfrog_matsuno') 
    130         CALL  leapfrog_matsuno_scheme 
    131  
    132       CASE ('adam_bashforth') 
    133         CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) 
    134         CALL adam_bashforth_scheme 
    135  
    136       CASE default 
    137         PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             & 
    138                ' > options are <euler>, <runge_kutta>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>' 
    139         STOP 
    140          
    141     END SELECT 
    142  
     159!    CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) 
    143160!    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q) 
    144      
    145  ENDDO 
    146    
     161    ENDDO 
     162  
    147163  CONTAINS 
    148164 
    … …  
    161177    END SUBROUTINE Euler_scheme 
    162178 
    163     SUBROUTINE RK_scheme 
     179    SUBROUTINE RK_scheme(stage) 
    164180      IMPLICIT NONE 
    165181      INTEGER :: ind, stage 
    166       REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ 1., 4./3., 2., 4. /) 
     182!      REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ 1., 4./3., 2., 4. /) 
     183      REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /) 
    167184      REAL(rstd) :: tau 
    168185 
    169       stage = MOD(it,4)+1 
    170186      tau = dt*coef(stage) 
    171  
     187       
    172188      DO ind=1,ndomain 
    173189         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind) 
    174190         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) 
    175191         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) 
    176  
    177          IF (stage==1) THEN 
     192          
     193         IF (stage==1) THEN ! first stage : save model state in XXm1 
    178194            psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) 
    179195         END IF 
     196         ! updates are of the form x1 := x0 + tau*f(x1) 
    180197         ps(:)=psm1(:)+tau*dps(:) 
    181198         u(:,:)=um1(:,:)+tau*du(:,:) 
    182199         theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 
    183200      END DO 
     201 
    184202    END SUBROUTINE RK_scheme 
    185203 
    … …  
    215233    END SUBROUTINE leapfrog_scheme   
    216234  
    217     SUBROUTINE leapfrog_matsuno_scheme 
     235    SUBROUTINE leapfrog_matsuno_scheme(stage) 
    218236    IMPLICIT NONE 
    219       INTEGER :: ind 
    220  
     237    INTEGER :: ind, stage 
     238    REAL :: tau 
     239    tau = dt/nb_stage 
    221240      DO ind=1,ndomain 
    222241        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind) 
    … …  
    226245 
    227246         
    228         IF (MOD(it,matsuno_period+1)==0) THEN 
     247!        IF (MOD(it,matsuno_period+1)==0) THEN 
     248        IF (stage==1) THEN 
    229249          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) 
    230           ps(:)=psm1(:)+dt*dps(:) 
    231           u(:,:)=um1(:,:)+dt*du(:,:) 
    232           theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:) 
    233  
    234         ELSE IF (MOD(it,matsuno_period+1)==1) THEN 
    235  
    236           ps(:)=psm1(:)+dt*dps(:) 
    237           u(:,:)=um1(:,:)+dt*du(:,:) 
    238           theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:) 
     250          ps(:)=psm1(:)+tau*dps(:) 
     251          u(:,:)=um1(:,:)+tau*du(:,:) 
     252          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 
     253 
     254!        ELSE IF (MOD(it,matsuno_period+1)==1) THEN 
     255        ELSE IF (stage==2) THEN 
     256 
     257          ps(:)=psm1(:)+tau*dps(:) 
     258          u(:,:)=um1(:,:)+tau*du(:,:) 
     259          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 
    239260 
    240261          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)  
    … …  
    243264        ELSE  
    244265 
    245           ps(:)=psm2(:)+2*dt*dps(:) 
    246           u(:,:)=um2(:,:)+2*dt*du(:,:) 
    247           theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:) 
     266          ps(:)=psm2(:)+2*tau*dps(:) 
     267          u(:,:)=um2(:,:)+2*tau*du(:,:) 
     268          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:) 
    248269 
    249270          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)  
Note: See TracChangeset for help on using the changeset viewer.