Changeset 360 for codes/icosagcm/trunk


Ignore:
Timestamp:
09/08/15 16:23:06 (9 years ago)
Author:
dubos
Message:

Towards HEVI time stepping

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

Legend:

Unmodified
Added
Removed

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

    r350 r360  
    5757  TYPE(t_field),POINTER :: f_hflux(:) 
    5858  TYPE(t_field),POINTER :: f_wflux(:) 
    59   TYPE(t_field),POINTER :: f_dps(:) 
    60   TYPE(t_field),POINTER :: f_dmass(:) 
    61   TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
    62   TYPE(t_field),POINTER :: f_du(:) 
     59  TYPE(t_field) :: f_dps(:) 
     60  TYPE(t_field) :: f_dmass(:) 
     61  TYPE(t_field) :: f_dtheta_rhodz(:) 
     62  TYPE(t_field) :: f_du(:) 
    6363 
    6464    SELECT CASE (TRIM(caldyn_type)) 
    … …  
    7171    END SELECT 
    7272   
    73   END SUBROUTINE caldyn   
     73  END SUBROUTINE caldyn 
    7474 
    7575END MODULE caldyn_mod 

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

    r343 r360  
    6464    TYPE(t_field),POINTER :: f_q(:) 
    6565    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:) 
    66     TYPE(t_field),POINTER :: f_dps(:) 
    67     TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
    68     TYPE(t_field),POINTER :: f_du(:) 
     66    TYPE(t_field) :: f_dps(:) 
     67    TYPE(t_field) :: f_dtheta_rhodz(:) 
     68    TYPE(t_field) :: f_du(:) 
    6969 
    7070    REAL(rstd),POINTER :: ps(:) 
    … …  
    103103  END SUBROUTINE caldyn 
    104104 
    105  
    106105  SUBROUTINE compute_caldyn(ps,u, hflux,wflux,dps, dtheta_rhodz,du) 
    107106    USE icosa 

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

    r354 r360  
    141141    TYPE(t_field),POINTER :: f_geopot(:) 
    142142    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:) 
    143     TYPE(t_field),POINTER :: f_dps(:) 
    144     TYPE(t_field),POINTER :: f_dmass(:) 
    145     TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
    146     TYPE(t_field),POINTER :: f_du(:) 
     143    TYPE(t_field) :: f_dps(:) 
     144    TYPE(t_field) :: f_dmass(:) 
     145    TYPE(t_field) :: f_dtheta_rhodz(:) 
     146    TYPE(t_field) :: f_du(:) 
    147147     
    148148    REAL(rstd),POINTER :: ps(:), dps(:) 

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

    r358 r360  
    241241    IF (is_master) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., & 
    242242                              'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000 
    243     IF (is_master)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=2.8 :', & 
    244                                2.8/340.*dumax**(-.5/nitergdiv) 
     243    IF (is_master)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=4 :', & 
     244                               4./340.*dumax**(-.5/nitergdiv) 
    245245 
    246246    cgraddiv=dumax**(-1./nitergdiv) 

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

    r354 r360  
    11MODULE timeloop_gcm_mod 
    2   USE transfert_mod 
    32  USE icosa 
     3  USE disvert_mod 
     4  USE trace 
     5  USE omp_para 
     6  USE euler_scheme_mod 
     7  USE explicit_scheme_mod 
     8  USE hevi_scheme_mod 
    49  IMPLICIT NONE 
    510  PRIVATE 
    611   
     12  INTEGER, PARAMETER :: itau_sync=10 
     13  TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0 
     14 
    715  PUBLIC :: init_timeloop, timeloop 
    8  
    9   INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3, rk25=4 
    10   INTEGER, PARAMETER :: itau_sync=10 
    11   REAL(rstd), DIMENSION(4), PARAMETER :: coef_rk4 = (/ .25, 1./3., .5, 1. /) 
    12   REAL(rstd), DIMENSION(5), PARAMETER :: coef_rk25 = (/ .25, 1./6., 3./8., .5, 1. /) 
    13  
    14   TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0 
    15  
    16   TYPE(t_field),POINTER,SAVE :: f_geopot(:) 
    17   TYPE(t_field),POINTER,SAVE :: f_q(:) 
    18   TYPE(t_field),POINTER,SAVE :: f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:) 
    19   TYPE(t_field),POINTER,SAVE :: f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:) 
    20   TYPE(t_field),POINTER,SAVE :: f_u(:),f_um1(:),f_um2(:), f_du(:) 
    21   TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:) 
    22   TYPE(t_field),POINTER,SAVE :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)   
    23  
    24   INTEGER,SAVE :: nb_stage, matsuno_period, scheme     
    25 !$OMP THREADPRIVATE(nb_stage, matsuno_period, scheme) 
    2616 
    2717CONTAINS 
    2818   
    2919  SUBROUTINE init_timeloop 
    30   USE icosa 
    31   USE dissip_gcm_mod 
    32   USE observable_mod 
    33   USE caldyn_mod 
    34   USE etat0_mod 
    35   USE disvert_mod 
    36   USE guided_mod 
    37   USE advect_tracer_mod 
    38   USE physics_mod 
    39   USE mpipara 
    40   USE omp_para 
    41   USE trace 
    42   USE transfert_mod 
    43   USE check_conserve_mod 
    44   USE output_field_mod 
    45   USE write_field_mod 
    46   USE theta2theta_rhodz_mod 
    47   USE sponge_mod 
    48  
    49   CHARACTER(len=255) :: def 
    50  
    51  
    52    IF (xios_output) itau_out=1 
    53    IF (.NOT. enable_io) itau_out=HUGE(itau_out) 
    54  
    55 ! Time-independant orography 
     20    USE dissip_gcm_mod 
     21    USE observable_mod 
     22    USE caldyn_mod 
     23    USE etat0_mod 
     24    USE guided_mod 
     25    USE advect_tracer_mod 
     26    USE check_conserve_mod 
     27    USE output_field_mod 
     28    USE theta2theta_rhodz_mod 
     29    USE sponge_mod 
     30 
     31    CHARACTER(len=255) :: def 
     32     
     33    IF (xios_output) itau_out=1 
     34    IF (.NOT. enable_io) itau_out=HUGE(itau_out) 
     35     
     36    def='RK2.5' 
     37    CALL getin('time_scheme',def) 
     38    SELECT CASE (TRIM(def)) 
     39    CASE('euler') 
     40       scheme_family=explicit 
     41       scheme=euler 
     42       nb_stage=1 
     43    CASE ('runge_kutta') 
     44       scheme_family=explicit 
     45       scheme=rk4 
     46       nb_stage=4 
     47    CASE ('RK2.5') 
     48       scheme_family=explicit 
     49       scheme=rk25 
     50       nb_stage=5 
     51    CASE ('leapfrog_matsuno') 
     52       scheme_family=explicit 
     53       scheme=mlf 
     54       matsuno_period=5 
     55       CALL getin('matsuno_period',matsuno_period) 
     56       nb_stage=matsuno_period+1 
     57    CASE('ARK2.3') 
     58       scheme_family=hevi 
     59       scheme=ark23 
     60       nb_stage=3 
     61       CALL set_coefs_ark23(dt) 
     62    CASE('ARK3.3') 
     63       scheme_family=hevi 
     64       scheme=ark33 
     65       nb_stage=3 
     66       CALL set_coefs_ark33(dt) 
     67    CASE ('none') 
     68       nb_stage=0 
     69    CASE default 
     70       PRINT*,'Bad selector for variable scheme : <', TRIM(def),             & 
     71            ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>,<RK2.5>,<ARK2.3>' 
     72       STOP 
     73    END SELECT 
     74     
     75    ! Time-independant orography 
    5676    CALL allocate_field(f_phis,field_t,type_real,name='phis') 
    57 ! Trends 
    58     CALL allocate_field(f_du,field_u,type_real,llm,name='du') 
    59     CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,name='dtheta_rhodz') 
    60 ! Model state at current time step (RK/MLF/Euler) 
     77    ! Model state at current time step 
     78    CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot') 
    6179    CALL allocate_field(f_ps,field_t,type_real, name='ps') 
    62     CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot') 
    6380    CALL allocate_field(f_mass,field_t,type_real,llm,name='mass') 
     81    CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz') 
     82    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz') 
    6483    CALL allocate_field(f_u,field_u,type_real,llm,name='u') 
    65     CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz') 
    66 ! Model state at previous time step (RK/MLF) 
    67     CALL allocate_field(f_um1,field_u,type_real,llm,name='um1') 
    68     CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,name='theta_rhodzm1') 
    69 ! Tracers 
    7084    CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q') 
    71     CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz') 
    72 ! Mass fluxes 
     85    ! Mass fluxes 
    7386    CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes 
    7487    CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time 
    7588    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! vertical mass fluxes 
    76     CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass') 
    77  
    78     IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default) 
     89    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt') 
     90     
     91    SELECT CASE(scheme_family) 
     92    CASE(explicit) 
     93       ! Trends 
    7994       CALL allocate_field(f_dps,field_t,type_real,name='dps') 
     95       CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass') 
     96       CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,name='dtheta_rhodz') 
     97       CALL allocate_field(f_du,field_u,type_real,llm,name='du') 
     98       ! Model state at previous time step (RK/MLF) 
    8099       CALL allocate_field(f_psm1,field_t,type_real,name='psm1') 
    81        CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt') 
    82        ! the following are unused but must point to something 
    83 !       f_massm1 => f_mass 
    84     ELSE ! eta = Lagrangian vertical coordinate 
    85100       CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1') 
    86        ! the following are unused but must point to something 
    87        f_wfluxt => f_wflux 
    88        f_dps  => f_phis 
    89        f_psm1 => f_phis 
    90     END IF 
    91  
    92     def='runge_kutta' 
    93     CALL getin('scheme',def) 
    94  
    95     SELECT CASE (TRIM(def)) 
    96       CASE('euler') 
    97         scheme=euler 
    98         nb_stage=1 
    99       CASE ('runge_kutta') 
    100         scheme=rk4 
    101         nb_stage=4 
    102       CASE ('RK2.5') 
    103         scheme=rk25 
    104         nb_stage=5 
    105       CASE ('leapfrog_matsuno') 
    106         scheme=mlf 
    107         matsuno_period=5 
    108         CALL getin('matsuno_period',matsuno_period) 
    109         nb_stage=matsuno_period+1 
    110      ! Model state 2 time steps ago (MLF) 
    111         CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm) 
    112         CALL allocate_field(f_um2,field_u,type_real,llm) 
    113         IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default) 
    114            CALL allocate_field(f_psm2,field_t,type_real) 
    115            ! the following are unused but must point to something 
    116            f_massm2 => f_mass 
    117         ELSE ! eta = Lagrangian vertical coordinate 
    118            CALL allocate_field(f_massm2,field_t,type_real,llm) 
    119            ! the following are unused but must point to something 
    120            f_psm2 => f_phis 
    121         END IF 
    122       CASE ('none') 
    123         nb_stage=0 
    124  
    125     CASE default 
    126        PRINT*,'Bad selector for variable scheme : <', TRIM(def),             & 
    127             ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>' 
    128        STOP 
     101       CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,name='theta_rhodzm1') 
     102       CALL allocate_field(f_um1,field_u,type_real,llm,name='um1') 
     103    CASE(hevi) 
     104       ! Trends 
     105       CALL allocate_fields(nb_stage,f_dps_slow, field_t,type_real,name='dps_slow') 
     106       CALL allocate_fields(nb_stage,f_dmass_slow, field_t,type_real,llm, name='dmass_slow') 
     107       CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,name='dtheta_rhodz_fast') 
     108       CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow') 
     109       CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast') 
     110       f_dps => f_dps_slow(:,1) 
     111       f_du => f_du_slow(:,1) 
     112       f_dtheta_rhodz => f_dtheta_rhodz_slow(:,1) 
     113    END SELECT 
     114 
     115    SELECT CASE(scheme) 
     116    CASE(mlf) 
     117       ! Model state 2 time steps ago (MLF) 
     118       CALL allocate_field(f_psm2,field_t,type_real) 
     119       CALL allocate_field(f_massm2,field_t,type_real,llm) 
     120       CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm) 
     121       CALL allocate_field(f_um2,field_u,type_real,llm) 
    129122    END SELECT 
    130123 
    … …  
    153146   
    154147  SUBROUTINE timeloop 
    155   USE icosa 
    156   USE dissip_gcm_mod 
    157   USE sponge_mod 
    158   USE disvert_mod 
    159   USE caldyn_mod 
    160   USE caldyn_gcm_mod, ONLY : req_ps, req_mass 
    161   USE observable_mod 
    162   USE etat0_mod 
    163   USE guided_mod 
    164   USE advect_tracer_mod 
    165   USE physics_mod 
    166   USE mpipara 
    167   USE omp_para 
    168   USE trace 
    169   USE transfert_mod 
    170   USE check_conserve_mod 
    171   USE xios_mod 
    172   USE output_field_mod 
    173   USE write_etat0_mod 
    174   USE checksum_mod 
    175   IMPLICIT NONE   
    176     REAL(rstd),POINTER :: q(:,:,:) 
    177     REAL(rstd),POINTER :: phis(:), ps(:) ,psm1(:), psm2(:), dps(:) 
    178     REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:), du(:,:) 
    179     REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), massm1(:,:), massm2(:,:), dmass(:,:) 
    180     REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:), dtheta_rhodz(:,:) 
    181     REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:) 
     148    USE dissip_gcm_mod 
     149    USE sponge_mod 
     150    USE observable_mod 
     151    USE etat0_mod 
     152    USE guided_mod 
     153    USE caldyn_mod 
     154    USE advect_tracer_mod 
     155    USE physics_mod 
     156    USE mpipara 
     157    USE transfert_mod 
     158    USE check_conserve_mod 
     159    USE xios_mod 
     160    USE output_field_mod 
     161    USE write_etat0_mod 
     162    USE checksum_mod 
     163    USE explicit_scheme_mod 
     164    USE hevi_scheme_mod 
     165    REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:) 
    182166 
    183167    INTEGER :: ind 
    184     INTEGER :: it,i,j,n,  stage 
     168    INTEGER :: it,i,j,l,n,  stage 
    185169    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 
    186     LOGICAL, PARAMETER :: check=.FALSE. 
    187     INTEGER :: start_clock 
    188     INTEGER :: stop_clock 
    189     INTEGER :: rate_clock 
    190     INTEGER :: l 
     170    LOGICAL, PARAMETER :: check_rhodz=.FALSE. 
     171    INTEGER :: start_clock, stop_clock, rate_clock 
    191172    LOGICAL,SAVE :: first_physic=.TRUE. 
    192 !$OMP THREADPRIVATE(first_physic)     
    193      
    194     
     173    !$OMP THREADPRIVATE(first_physic)     
     174 
    195175    CALL switch_omp_distrib_level 
    196176    CALL caldyn_BC(f_phis, f_geopot, f_wflux) ! set constant values in first/last interfaces 
    197      
    198 !$OMP BARRIER 
    199   DO ind=1,ndomain 
    200      IF (.NOT. assigned_domain(ind)) CYCLE 
    201      CALL swap_dimensions(ind) 
    202      CALL swap_geometry(ind) 
    203      rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind) 
    204      IF(caldyn_eta==eta_mass) THEN 
    205         CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps 
    206      ELSE 
    207        DO l=ll_begin,ll_end 
    208          rhodz(:,l)=mass(:,l) 
    209        ENDDO 
    210      END IF 
    211   END DO 
    212 !$OMP BARRIER 
    213   fluxt_zero=.TRUE. 
    214  
    215 !$OMP MASTER 
    216   CALL SYSTEM_CLOCK(start_clock) 
    217 !$OMP END MASTER    
    218  
    219   CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)   
    220  
    221   CALL trace_on 
    222    
    223   DO it=itau0+1,itau0+itaumax 
    224       
    225     CALL check_conserve_detailed('detailed_budget 0', & 
    226          f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
    227  
    228     IF (xios_output) CALL xios_update_calendar(it) 
    229     IF (it==itau0+1 .OR. MOD(it,itau_sync)==0) THEN 
    230       CALL send_message(f_ps,req_ps0) 
    231       CALL wait_message(req_ps0) 
    232       CALL send_message(f_mass,req_mass0) 
    233       CALL wait_message(req_mass0) 
    234       CALL send_message(f_theta_rhodz,req_theta_rhodz0)  
    235       CALL wait_message(req_theta_rhodz0)  
    236       CALL send_message(f_u,req_u0) 
    237       CALL wait_message(req_u0) 
    238       CALL send_message(f_q,req_q0)  
    239       CALL wait_message(req_q0)  
    240  
    241     ENDIF 
    242  
    243     IF (is_master) PRINT *,"It No :",It,"   t :",dt*It 
    244  
    245     IF (mod(it,itau_out)==0 ) THEN 
    246       CALL update_time_counter(dt*it) 
    247       CALL write_output_fields_basic(f_ps, f_u, f_q) 
    248 !      CALL output_field("q",f_q) 
    249     ENDIF 
    250      
    251     CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 
    252  
    253     DO stage=1,nb_stage 
    254 !       CALL checksum(f_ps) 
    255 !       CALL checksum(f_theta_rhodz) 
    256 !       CALL checksum(f_mass) 
    257        CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), & 
    258             f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, & 
    259             f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du) 
    260 !       CALL checksum(f_dps) 
    261 !       CALL checksum(f_dtheta_rhodz) 
    262 !       CALL checksum(f_dmass) 
    263        SELECT CASE (scheme) 
    264        CASE(euler) 
    265           CALL euler_scheme(.TRUE.) 
    266        CASE (rk4) 
    267           CALL rk_scheme(stage, coef_rk4) 
    268        CASE (rk25) 
    269           CALL rk_scheme(stage, coef_rk25) 
    270        CASE (mlf) 
    271           CALL  leapfrog_matsuno_scheme(stage) 
    272        CASE DEFAULT 
    273           STOP 
    274        END SELECT 
    275     END DO 
    276  
    277     CALL check_conserve_detailed('detailed_budget 1', & 
    278          f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
    279  
    280     IF (MOD(it,itau_dissip)==0) THEN 
    281         
    282        IF(caldyn_eta==eta_mass) THEN 
    283 !ym flush ps 
    284 !$OMP BARRIER 
    285           DO ind=1,ndomain 
    286              IF (.NOT. assigned_domain(ind)) CYCLE 
    287              CALL swap_dimensions(ind) 
    288              CALL swap_geometry(ind) 
    289              mass=f_mass(ind); ps=f_ps(ind); 
    290              CALL compute_rhodz(.TRUE., ps, mass) 
    291           END DO 
    292        ENDIF 
    293  
    294        CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz) 
    295  
    296        CALL euler_scheme(.FALSE.)  ! update only u, theta 
    297        IF (iflag_sponge > 0) THEN 
    298          CALL sponge(f_u,f_du,f_theta_rhodz,f_dtheta_rhodz) 
    299          CALL euler_scheme(.FALSE.)  ! update only u, theta 
    300        ENDIF 
    301     END IF 
    302  
    303     CALL check_conserve_detailed('detailed_budget 2', & 
    304          f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
    305  
    306     IF(MOD(it,itau_adv)==0) THEN 
    307  
    308        CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step 
    309        fluxt_zero=.TRUE. 
    310  
    311        ! FIXME : check that rhodz is consistent with ps 
    312        IF (check) THEN 
    313          DO ind=1,ndomain 
    314             IF (.NOT. assigned_domain(ind)) CYCLE 
    315             CALL swap_dimensions(ind) 
    316             CALL swap_geometry(ind) 
    317             rhodz=f_rhodz(ind); ps=f_ps(ind); 
    318             CALL compute_rhodz(.FALSE., ps, rhodz)    
    319          END DO 
    320        ENDIF 
    321  
    322     END IF 
    323  
    324  
    325     IF (MOD(it,itau_check_conserv)==0) THEN 
    326       CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
    327     ENDIF 
    328       
    329     IF (MOD(it,itau_physics)==0) THEN 
    330       CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q) 
    331  
    332 !$OMP MASTER 
    333    IF (first_physic) CALL SYSTEM_CLOCK(start_clock) 
    334 !$OMP END MASTER    
    335       first_physic=.FALSE. 
    336     ENDIF 
    337      
    338   ENDDO 
    339  
    340   CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q)  
    341  
    342   CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)   
    343  
    344 !$OMP MASTER 
    345   CALL SYSTEM_CLOCK(stop_clock) 
    346   CALL SYSTEM_CLOCK(count_rate=rate_clock) 
    347      
    348   IF (mpi_rank==0) THEN  
    349     PRINT *,"Time elapsed : ",(stop_clock-start_clock)*1./rate_clock  
    350   ENDIF   
    351 !$OMP END MASTER  
    352   
    353  CONTAINS 
    354  
    355     SUBROUTINE Euler_scheme(with_dps) 
    356     IMPLICIT NONE 
    357     LOGICAL :: with_dps 
    358     INTEGER :: ind 
    359     INTEGER :: i,j,ij,l 
    360     CALL trace_start("Euler_scheme")   
    361  
     177 
     178    !$OMP BARRIER 
    362179    DO ind=1,ndomain 
    363180       IF (.NOT. assigned_domain(ind)) CYCLE 
    364181       CALL swap_dimensions(ind) 
    365182       CALL swap_geometry(ind) 
    366  
    367        IF(with_dps) THEN ! update ps/mass 
    368           IF(caldyn_eta==eta_mass) THEN ! update ps 
    369              ps=f_ps(ind) ; dps=f_dps(ind) ;               
    370              IF (is_omp_first_level) THEN 
    371 !$SIMD 
    372                 DO ij=ij_begin,ij_end 
    373                     ps(ij)=ps(ij)+dt*dps(ij) 
    374                 ENDDO 
    375              ENDIF  
    376           ELSE ! update mass 
    377              mass=f_mass(ind) ; dmass=f_dmass(ind) ;               
    378              DO l=ll_begin,ll_end 
    379 !$SIMD 
    380                 DO ij=ij_begin,ij_end 
    381                     mass(ij,l)=mass(ij,l)+dt*dmass(ij,l) 
    382                 ENDDO 
    383              END DO 
    384           END IF 
    385  
    386           hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind) 
    387           wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind) 
    388           CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind)) 
    389        END IF ! update ps/mass 
    390         
    391        u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) 
    392        du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) 
    393  
    394        DO l=ll_begin,ll_end 
    395 !$SIMD 
    396          DO ij=ij_begin,ij_end 
    397              u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l) 
    398              u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l) 
    399              u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l) 
    400              theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l) 
    401          ENDDO 
    402        ENDDO 
    403     ENDDO 
    404  
    405     CALL trace_end("Euler_scheme")   
    406  
    407     END SUBROUTINE Euler_scheme 
    408  
    409     SUBROUTINE RK_scheme(stage,coef) 
    410       IMPLICIT NONE 
    411       INTEGER :: ind, stage 
    412       REAL(rstd), INTENT(IN) :: coef(:) 
    413       REAL(rstd) :: tau 
    414       INTEGER :: i,j,ij,l 
    415    
    416       CALL trace_start("RK_scheme")   
    417  
    418       tau = dt*coef(stage) 
    419  
    420       ! if mass coordinate, deal with ps first on one core 
    421       IF(caldyn_eta==eta_mass) THEN 
    422          IF (is_omp_first_level) THEN 
    423  
    424             DO ind=1,ndomain 
    425                IF (.NOT. assigned_domain(ind)) CYCLE 
    426                CALL swap_dimensions(ind) 
    427                CALL swap_geometry(ind) 
    428                ps=f_ps(ind) ; psm1=f_psm1(ind) ; dps=f_dps(ind)  
    429                 
    430                IF (stage==1) THEN ! first stage : save model state in XXm1 
    431 !$SIMD 
    432                  DO ij=ij_begin,ij_end 
    433                    psm1(ij)=ps(ij) 
    434                  ENDDO 
    435                ENDIF 
    436              
    437                ! updates are of the form x1 := x0 + tau*f(x1) 
    438 !$SIMD 
    439                DO ij=ij_begin,ij_end 
    440                    ps(ij)=psm1(ij)+tau*dps(ij) 
    441                ENDDO 
    442             ENDDO 
    443          ENDIF 
    444 !         CALL send_message(f_ps,req_ps) 
    445 !ym no overlap for now          
    446 !         CALL wait_message(req_ps)   
    447        
    448       ELSE ! Lagrangian coordinate, deal with mass 
    449          DO ind=1,ndomain 
    450             IF (.NOT. assigned_domain(ind)) CYCLE 
    451             CALL swap_dimensions(ind) 
    452             CALL swap_geometry(ind) 
    453             mass=f_mass(ind); dmass=f_dmass(ind); massm1=f_massm1(ind) 
    454  
    455             IF (stage==1) THEN ! first stage : save model state in XXm1 
    456                DO l=ll_begin,ll_end 
    457 !$SIMD 
    458                  DO ij=ij_begin,ij_end 
    459                     massm1(ij,l)=mass(ij,l) 
    460                  ENDDO 
    461                ENDDO 
    462             END IF 
    463  
    464             ! updates are of the form x1 := x0 + tau*f(x1) 
    465             DO l=ll_begin,ll_end 
    466 !$SIMD 
    467                DO ij=ij_begin,ij_end 
    468                    mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l) 
    469                ENDDO 
    470             ENDDO 
    471          END DO 
    472 !         CALL send_message(f_mass,req_mass) 
    473 !ym no overlap for now          
    474 !         CALL wait_message(req_mass)   
    475  
    476       END IF 
    477  
    478       ! now deal with other prognostic variables 
    479       DO ind=1,ndomain 
    480          IF (.NOT. assigned_domain(ind)) CYCLE 
    481          CALL swap_dimensions(ind) 
    482          CALL swap_geometry(ind) 
    483          u=f_u(ind)      ; du=f_du(ind)      ; um1=f_um1(ind)  
    484          theta_rhodz=f_theta_rhodz(ind) 
    485          theta_rhodzm1=f_theta_rhodzm1(ind) 
    486          dtheta_rhodz=f_dtheta_rhodz(ind) 
    487           
    488          IF (stage==1) THEN ! first stage : save model state in XXm1 
    489            DO l=ll_begin,ll_end 
    490 !$SIMD 
    491                 DO ij=ij_begin,ij_end 
    492                  um1(ij+u_right,l)=u(ij+u_right,l) 
    493                  um1(ij+u_lup,l)=u(ij+u_lup,l) 
    494                  um1(ij+u_ldown,l)=u(ij+u_ldown,l) 
    495                  theta_rhodzm1(ij,l)=theta_rhodz(ij,l) 
    496              ENDDO 
    497            ENDDO 
    498          END IF         
    499  
    500          DO l=ll_begin,ll_end 
    501 !$SIMD 
    502            DO ij=ij_begin,ij_end 
    503                u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l) 
    504                u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l) 
    505                u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l) 
    506                theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l) 
    507            ENDDO 
    508          ENDDO 
    509  
    510          IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage 
    511             hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind) 
    512             wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind) 
    513             CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind)) 
    514          END IF 
    515       END DO 
    516        
    517       CALL trace_end("RK_scheme") 
    518        
    519     END SUBROUTINE RK_scheme 
    520  
    521     SUBROUTINE leapfrog_matsuno_scheme(stage) 
    522     IMPLICIT NONE 
    523     INTEGER :: ind, stage 
    524     REAL :: tau 
    525  
    526       CALL trace_start("leapfrog_matsuno_scheme") 
    527      
    528       tau = dt/nb_stage 
    529       DO ind=1,ndomain 
    530         IF (.NOT. assigned_domain(ind)) CYCLE 
    531         CALL swap_dimensions(ind) 
    532         CALL swap_geometry(ind) 
    533  
    534         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind) 
    535         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) 
    536         psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind) 
    537         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) 
    538  
    539          
    540         IF (stage==1) THEN 
    541           psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) 
    542           ps(:)=psm1(:)+tau*dps(:) 
    543           u(:,:)=um1(:,:)+tau*du(:,:) 
    544           theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 
    545  
    546         ELSE IF (stage==2) THEN 
    547  
    548           ps(:)=psm1(:)+tau*dps(:) 
    549           u(:,:)=um1(:,:)+tau*du(:,:) 
    550           theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 
    551  
    552           psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)  
    553           psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)  
    554  
    555         ELSE  
    556  
    557           ps(:)=psm2(:)+2*tau*dps(:) 
    558           u(:,:)=um2(:,:)+2*tau*du(:,:) 
    559           theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:) 
    560  
    561           psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)  
    562           psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)  
    563  
    564         ENDIF 
    565        
    566       ENDDO 
    567       CALL trace_end("leapfrog_matsuno_scheme") 
    568        
    569     END SUBROUTINE leapfrog_matsuno_scheme   
    570           
    571   END SUBROUTINE timeloop     
    572  
    573  SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero) 
    574     USE icosa 
    575     USE omp_para 
    576     USE disvert_mod 
    577     IMPLICIT NONE 
    578     REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1) 
    579     REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1) 
    580     REAL(rstd), INTENT(IN) :: tau 
    581     LOGICAL, INTENT(INOUT) :: fluxt_zero 
    582     INTEGER :: l,i,j,ij 
    583  
    584     IF(fluxt_zero) THEN 
    585  
    586        fluxt_zero=.FALSE. 
    587  
    588        DO l=ll_begin,ll_end 
    589 !$SIMD 
    590          DO ij=ij_begin_ext,ij_end_ext 
    591              hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l) 
    592              hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l) 
    593              hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l) 
    594          ENDDO 
    595        ENDDO 
    596  
    597        IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 
    598           DO l=ll_begin,ll_endp1 
    599 !$SIMD 
    600              DO ij=ij_begin,ij_end 
    601                 wfluxt(ij,l) = tau*wflux(ij,l) 
    602              ENDDO 
     183       rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind) 
     184       IF(caldyn_eta==eta_mass) THEN 
     185          CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps 
     186       ELSE 
     187          DO l=ll_begin,ll_end 
     188             rhodz(:,l)=mass(:,l) 
    603189          ENDDO 
    604190       END IF 
    605  
    606     ELSE 
    607  
    608        DO l=ll_begin,ll_end 
    609 !$SIMD 
    610          DO ij=ij_begin_ext,ij_end_ext 
    611              hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l) 
    612              hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l) 
    613              hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l) 
    614          ENDDO 
    615        ENDDO 
    616  
    617        IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 
    618           DO l=ll_begin,ll_endp1 
    619 !$SIMD 
    620              DO ij=ij_begin,ij_end 
    621                    wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l) 
    622              ENDDO 
    623           ENDDO 
     191    END DO 
     192    !$OMP BARRIER 
     193    fluxt_zero=.TRUE. 
     194 
     195    !$OMP MASTER 
     196    CALL SYSTEM_CLOCK(start_clock) 
     197    !$OMP END MASTER    
     198 
     199    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)   
     200 
     201    CALL trace_on 
     202 
     203    DO it=itau0+1,itau0+itaumax 
     204 
     205       CALL check_conserve_detailed('detailed_budget 0', & 
     206            f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
     207 
     208       IF (xios_output) CALL xios_update_calendar(it) 
     209       IF (it==itau0+1 .OR. MOD(it,itau_sync)==0) THEN 
     210          CALL send_message(f_ps,req_ps0) 
     211          CALL wait_message(req_ps0) 
     212          CALL send_message(f_mass,req_mass0) 
     213          CALL wait_message(req_mass0) 
     214          CALL send_message(f_theta_rhodz,req_theta_rhodz0)  
     215          CALL wait_message(req_theta_rhodz0)  
     216          CALL send_message(f_u,req_u0) 
     217          CALL wait_message(req_u0) 
     218          CALL send_message(f_q,req_q0)  
     219          CALL wait_message(req_q0)  
     220       ENDIF 
     221 
     222       IF (is_master) PRINT *,"It No :",It,"   t :",dt*It 
     223 
     224       IF (mod(it,itau_out)==0 ) THEN 
     225          CALL update_time_counter(dt*it) 
     226          CALL write_output_fields_basic(f_ps, f_u, f_q) 
     227       ENDIF 
     228 
     229       CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 
     230 
     231       SELECT CASE(scheme_family) 
     232       CASE(explicit) 
     233          CALL explicit_scheme(it, fluxt_zero) 
     234       CASE(hevi) 
     235          CALL HEVI_scheme(it, fluxt_zero) 
     236       END SELECT 
     237 
     238       CALL check_conserve_detailed('detailed_budget 1', & 
     239            f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
     240        
     241       IF (MOD(it,itau_dissip)==0) THEN 
     242           
     243          IF(caldyn_eta==eta_mass) THEN 
     244             !ym flush ps 
     245             !$OMP BARRIER 
     246             DO ind=1,ndomain 
     247                IF (.NOT. assigned_domain(ind)) CYCLE 
     248                CALL swap_dimensions(ind) 
     249                CALL swap_geometry(ind) 
     250                mass=f_mass(ind); ps=f_ps(ind); 
     251                CALL compute_rhodz(.TRUE., ps, mass) 
     252             END DO 
     253          ENDIF 
     254           
     255          CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz) 
     256           
     257          CALL euler_scheme(.FALSE.)  ! update only u, theta 
     258          IF (iflag_sponge > 0) THEN 
     259             CALL sponge(f_u,f_du,f_theta_rhodz,f_dtheta_rhodz) 
     260             CALL euler_scheme(.FALSE.)  ! update only u, theta 
     261          END IF 
    624262       END IF 
    625  
    626    END IF 
    627  
    628   END SUBROUTINE accumulate_fluxes 
     263        
     264       CALL check_conserve_detailed('detailed_budget 2', & 
     265            f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
     266        
     267       IF(MOD(it,itau_adv)==0) THEN 
     268          CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step 
     269          fluxt_zero=.TRUE. 
     270          ! FIXME : check that rhodz is consistent with ps 
     271          IF (check_rhodz) THEN 
     272             DO ind=1,ndomain 
     273                IF (.NOT. assigned_domain(ind)) CYCLE 
     274                CALL swap_dimensions(ind) 
     275                CALL swap_geometry(ind) 
     276                rhodz=f_rhodz(ind); ps=f_ps(ind); 
     277                CALL compute_rhodz(.FALSE., ps, rhodz)    
     278             END DO 
     279          ENDIF 
     280           
     281       END IF 
     282        
     283       IF (MOD(it,itau_check_conserv)==0) THEN 
     284          CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
     285       ENDIF 
     286        
     287       IF (MOD(it,itau_physics)==0) THEN 
     288          CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)         
     289          !$OMP MASTER 
     290          IF (first_physic) CALL SYSTEM_CLOCK(start_clock) 
     291          !$OMP END MASTER    
     292          first_physic=.FALSE. 
     293       END IF 
     294        
     295    END DO 
     296     
     297    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q)  
     298     
     299    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)   
     300     
     301    !$OMP MASTER 
     302    CALL SYSTEM_CLOCK(stop_clock) 
     303    CALL SYSTEM_CLOCK(count_rate=rate_clock) 
     304     
     305    IF (mpi_rank==0) THEN  
     306       PRINT *,"Time elapsed : ",(stop_clock-start_clock)*1./rate_clock  
     307    ENDIF 
     308    !$OMP END MASTER  
     309     
     310    ! CONTAINS 
     311  END SUBROUTINE timeloop 
    629312   
    630313END MODULE timeloop_gcm_mod 
Note: See TracChangeset for help on using the changeset viewer.