Changeset 642


Ignore:
Timestamp:
12/19/17 15:26:51 (6 years ago)
Author:
dubos
Message:

devel/unstructured : bubble test case with Fortran time stepping

Location:
codes/icosagcm/devel
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/devel/Python/dynamico/meshes.py

    r639 r642  
    148148                  Aiv,Aiv,f+0.*Aiv,le_de,Riv2,-wee) 
    149149 
    150     def field_theta(self): return squeeze((self.nqdyn,self.ny,self.nx,self.llm)) 
    151     def field_mass(self): return squeeze((self.ny,self.nx,self.llm)) 
    152     def field_z(self): return squeeze((self.ny,self.nx,self.llm)) 
    153     def field_w(self): return squeeze((self.ny,self.nx,self.llm+1)) 
    154     def field_u(self): return np.zeros((self.ny,2*self.nx,self.llm)) 
    155     def field_ps(self): return squeeze((self.ny,self.nx)) 
     150    def field_theta(self,n=1): return squeeze((n,self.nqdyn,self.ny,self.nx,self.llm)) 
     151    def field_mass(self,n=1): return squeeze((n,self.ny,self.nx,self.llm)) 
     152    def field_z(self,n=1): return squeeze((n,self.ny,self.nx,self.llm)) 
     153    def field_w(self,n=1): return squeeze((n,self.ny,self.nx,self.llm+1)) 
     154    def field_u(self,n=1):  
     155        if n==1: 
     156            return np.zeros((self.ny,2*self.nx,self.llm)) 
     157        else: 
     158            return np.zeros((n,self.ny,2*self.nx,self.llm)) 
     159    def field_ps(self,n=1): return squeeze((n,self.ny,self.nx)) 
    156160    def ucomp(self,u): return u[:,range(0,2*self.nx,2),:] 
    157161    def set_ucomp(self,uv,u): uv[:,range(0,2*self.nx,2),:]=u 

  • codes/icosagcm/devel/Python/src/functions.h

    r639 r642  
    77void dynamico_init_params(void); 
    88 
    9 void dynamico_ARK_step(double *mass_col, double *rhodz, double *theta_rhodz,  
     9void dynamico_ARK_step(int nstep, 
     10                       double *mass_col, double *rhodz, double *theta_rhodz,  
    1011                       double *u, double *geopot, double *w, 
    1112                       double *theta, double *ps, double *pk, double *hflux, double *qv, 

  • codes/icosagcm/devel/Python/src/unstructured.pyx

    r640 r642  
    1818 
    1919cdef extern from "functions.h": 
    20      cdef void dynamico_ARK_step(double *mass_col, double *rhodz, double *theta_rhodz,  
     20     cdef void dynamico_ARK_step(int nstep, 
     21                                 double *mass_col, double *rhodz, double *theta_rhodz,  
    2122                                 double *u, double *geopot, double *w, 
    2223                                 double *theta, double *ps, double *pk, double *hflux, double *qv, 
    … …  
    7879cdef dbl_ptr ptr2(double[:,:] data) : return &data[0,0] 
    7980cdef dbl_ptr ptr3(double[:,:,:] data) : return &data[0,0,0] 
     81cdef dbl_ptr ptr4(double[:,:,:,:] data) : return &data[0,0,0,0] 
    8082cdef dbl_ptr ptr(data): 
    8183    n=data.ndim 
    … …  
    8385    if n==2 : return ptr2(data) 
    8486    if n==3 : return ptr3(data) 
    85  
     87    if n==4 : return ptr4(data) 
     88    if n>4: raise IndexError 
     89         
    8690cdef alloc(dbl_ptr *p, allocator, n=1): 
    8791    data=allocator(n) 
    … …  
    9397         
    9498cdef class Caldyn_step: 
     99    # number of time steps to do at each invocation of advance() 
     100    cdef int nstep 
    95101    # pointer to allocated arrays 
    96     cdef dbl_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_w                   # prognostic 
     102    cdef dbl_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_W                   # prognostic 
    97103    cdef dbl_ptr p_mass_col, p_dmass_col, p_ps, p_theta, p_pk, p_hflux, p_qv # diagnostic 
    98104    cdef dbl_ptr p_drhodz, p_dtheta_rhodz, p_du_fast, p_du_slow                # tendencies 
    99105    cdef dbl_ptr p_dPhi_fast, p_dPhi_slow, p_dW_fast, p_dW_slow                # tendencies 
    100106    # allocated arrays, must remain referenced or segfault 
    101     cdef readonly np.ndarray mass, theta_rhodz, u, geopot, w 
     107    cdef readonly np.ndarray mass, theta_rhodz, u, geopot, W 
    102108    cdef readonly np.ndarray mass_col, dmass_col, ps, theta, pk, hflux, qv 
    103109    cdef readonly np.ndarray drhodz, dtheta_rhodz, du_fast, du_slow 
    104110    cdef readonly np.ndarray dPhi_fast, dPhi_slow, dW_fast, dW_slow 
    105111 
    106     def __init__(self,mesh,time_scheme): 
     112    def __init__(self,mesh,time_scheme, nstep): 
     113        self.nstep=nstep 
    107114        #        self.mesh=mesh 
    108115        fps, ftheta, fmass = mesh.field_ps, mesh.field_theta, mesh.field_mass 
    … …  
    113120        cdef double[:,:] cflj_ = time_scheme.cfjl 
    114121        ns = time_scheme.nstage 
    115         nb_stage[0]= ns 
     122        nb_stage[0]=ns 
     123 
    116124        cdef int i,j 
    117125        for i in range(ns): 
    … …  
    123131        #    prognostic/diagnostic 
    124132        self.ps                       = alloc(&self.p_ps, fps)  
    125         self.mass_col, self.dmass_col = alloc(&self.p_mass_col, fps), alloc(&self.p_dmass_col, fps),  
     133        self.mass_col, self.dmass_col = alloc(&self.p_mass_col, fps), alloc(&self.p_dmass_col, fps,ns),  
    126134        self.mass, self.theta_rhodz   = alloc(&self.p_mass, fmass), alloc(&self.p_theta_rhodz, fmass), 
    127135        self.theta, self.pk           = alloc(&self.p_theta, fmass), alloc(&self.p_pk, fmass),  
    128         self.geopot, self.w           = alloc(&self.p_geopot, fw), alloc(&self.p_w, fw), 
     136        self.geopot, self.W           = alloc(&self.p_geopot, fw), alloc(&self.p_W, fw), 
    129137        self.hflux, self.u            = alloc(&self.p_hflux, fu), alloc(&self.p_u, fu)  
    130138        self.qv                       = alloc(&self.p_qv,fz) 
    … …  
    137145        #        global elapsed 
    138146        # time1=time.time() 
    139         check_ptr('mass', self.p_mass, self.mass) 
    140         check_ptr('u', self.p_u, self.u) 
    141         dynamico_ARK_step(self.p_mass_col, self.p_mass, self.p_theta_rhodz, 
    142                           self.p_u, self.p_geopot, self.p_w, 
     147        dynamico_ARK_step(self.nstep, 
     148                          self.p_mass_col, self.p_mass, self.p_theta_rhodz, 
     149                          self.p_u, self.p_geopot, self.p_W, 
    143150                          self.p_theta, self.p_ps, self.p_pk, self.p_hflux, self.p_qv, 
    144151                          self.p_dmass_col, self.p_drhodz, self.p_dtheta_rhodz, 
    … …  
    148155        #time2=time.time() 
    149156        #if time2>time1: elapsed=elapsed+time2-time1 
    150     def setup_TRSW(self): 
    151         setvars(('hydrostatic','caldyn_thermo','caldyn_eta'), 
    152                 (True,thermo_boussinesq,eta_lag)) 
    153         dynamico_init_params() 
    154     def setup_HPE(self, caldyn_thermo, caldyn_eta, g, BC,thermo): 
    155         setvars(('hydrostatic','caldyn_thermo','caldyn_eta', 
    156                  'g','ptop','Rd','cpp','preff','Treff'), 
    157                 (True,caldyn_thermo,caldyn_eta, 
    158                  g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0)) 
    159         dynamico_init_params() 
    160     def setup_NH(self, caldyn_thermo, caldyn_eta, g, BC,thermo): 
    161         setvars(('hydrostatic','caldyn_thermo','caldyn_eta', 
    162                  'g','ptop','Rd','cpp','preff','Treff', 
    163                  'pbot','rho_bot'), 
    164                 (False,caldyn_thermo,caldyn_eta, 
    165                  g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0, 
    166                  BC.pbot.max(), BC.rho_bot.max())) 
    167         dynamico_init_params() 
    168  
     157 
     158def caldyn_step_TRSW(mesh,time_scheme,nstep): 
     159    setvars(('hydrostatic','caldyn_thermo','caldyn_eta'), 
     160            (True,thermo_boussinesq,eta_lag)) 
     161    dynamico_init_params() 
     162    return Caldyn_step(mesh,time_scheme, nstep) 
     163def caldyn_step_HPE(mesh,time_scheme,nstep, caldyn_thermo,caldyn_eta, thermo,BC,g): 
     164    setvars(('hydrostatic','caldyn_thermo','caldyn_eta', 
     165             'g','ptop','Rd','cpp','preff','Treff'), 
     166             (True,caldyn_thermo,caldyn_eta, 
     167              g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0)) 
     168    dynamico_init_params() 
     169    return Caldyn_step(mesh,time_scheme, nstep) 
     170def caldyn_step_NH(mesh,time_scheme,nstep, caldyn_thermo, caldyn_eta, thermo,BC,g): 
     171    setvars(('hydrostatic','caldyn_thermo','caldyn_eta', 
     172             'g','ptop','Rd','cpp','preff','Treff','pbot','rho_bot'), 
     173             (False,caldyn_thermo,caldyn_eta, 
     174              g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0, 
     175              BC.pbot.max(), BC.rho_bot.max())) 
     176    dynamico_init_params() 
     177    return Caldyn_step(mesh,time_scheme, nstep) 
     178     
    169179#----------------------------- Base class for dynamics ------------------------ 
    170180 
    … …  
    212222 
    213223class Caldyn_HPE(Caldyn): 
    214     def __init__(self,caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g): 
     224    def __init__(self,caldyn_thermo,caldyn_eta, mesh,thermo,BC,g): 
    215225        Caldyn.__init__(self,mesh) 
    216226        setvars(('hydrostatic','caldyn_thermo','caldyn_eta', 
    … …  
    231241 
    232242class Caldyn_NH(Caldyn): 
    233     def __init__(self,caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g): 
     243    def __init__(self,caldyn_thermo,caldyn_eta, mesh,thermo,BC,g): 
    234244        Caldyn.__init__(self,mesh) 
    235245        setvars(('hydrostatic','caldyn_thermo','caldyn_eta', 

  • codes/icosagcm/devel/Python/test/fig_RSW2_MPAS_W02

    • Property svn:ignore set to
      *.png

  • codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py

    r632 r642  
    5353    return thermo, mesh, params, (mass_ik,Sik,ujk,Phi_il,Wil), gas 
    5454 
    55 Lx, nx, llm, thetac, T, Nslice, courant = 2000., 100, 50, 30., 5., 10, 2.8 
    56 #Lx, nx, llm, thetac, T, Nslice, courant = 2000., 50, 25, 1., 50., 10, 2.8 
     55#Lx, nx, llm, thetac, T, Nslice, courant = 2000., 100, 50, 30., 5., 10, 2.8 
     56Lx, nx, llm, thetac, T, Nslice, courant = 2000., 50, 25, 30, 5., 10, 2.8 
    5757#Lx, nx, llm, thetac, T, Nslice, courant = 3000., 75, 25, -30, 5., 10, 2.8 
    5858 
    … …  
    6161thermo, mesh, params, flow0, gas0 = thermal_bubble_3D(Lx,nx,Ly,ny,llm, thetac=thetac) 
    6262 
    63 caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass 
    64 caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,params,thermo,params,params.g) 
    65  
    6663# compute hybrid coefs from initial distribution of mass 
    6764mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0]) 
    6865unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk) 
    69 print mass_dbk 
    7066 
    7167dz = flow0[3].max()/(params.g*llm) 
    … …  
    7470nt = int(math.ceil(T/dt)) 
    7571dt = T/nt 
    76 print 'Time step : %g s' % dt 
     72print 'Time step : %d x %g s' % (nt,dt) 
    7773 
    78 scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt) 
     74caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass 
     75#caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_lag 
    7976 
    80 flow=flow0 
     77if False: # time stepping in Python 
     78    caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,thermo,params,params.g) 
     79    scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt) 
     80    def next_flow(m,S,u,Phi,W): 
     81        # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) 
     82        return scheme.advance((m,S,u,Phi,W),nt) 
     83 
     84else: # time stepping in Fortran 
     85    scheme = time_step.ARK2(None, dt) 
     86    caldyn_step = unst.caldyn_step_NH(mesh,scheme,nt, caldyn_thermo,caldyn_eta, 
     87                                      thermo,params,params.g) 
     88    def next_flow(m,S,u,Phi,W): 
     89        # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) 
     90        caldyn_step.mass[:,:,:], caldyn_step.theta_rhodz[:,:,:], caldyn_step.u[:,:,:] = m,S,u 
     91        caldyn_step.geopot[:,:,:], caldyn_step.W[:,:,:] = Phi,W 
     92        caldyn_step.next() 
     93        return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(), 
     94                caldyn_step.geopot.copy(), caldyn_step.W.copy()) 
     95     
     96m,S,u,Phi,W=flow0 
    8197w=mesh.field_mass() 
    8298z=mesh.field_mass() 
     99 
    83100for it in range(Nslice): 
    84     junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) 
    85     m,S,u,Phi,W = flow 
    86     s=S/m ; # s=.5*(s+abs(s)) 
     101    s=S/m ; s=.5*(s+abs(s)) 
    87102    for l in range(llm): 
    88103        w[:,:,l]=.5*params.g*(W[:,:,l+1]+W[:,:,l])/m[:,:,l] 
    … …  
    92107    jj=ny/2 
    93108    xx,zz,ss,ww = mesh.xx[:,jj,:]/1000., z[:,jj,:]/1000., s[:,jj,:], w[:,jj,:] 
     109 
     110    f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12,4)) 
    94111     
    95     f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12,4)) 
    96  
    97112    c=ax1.contourf(xx,zz,ss,20) 
    98113    ax1.set_xlim((-.5,.5)), ax1.set_xlabel('x (km)') 
    … …  
    100115    plt.colorbar(c,ax=ax1) 
    101116    ax1.set_title('Specific entropy at t=%g s (J/K/kg)' % (it*T,)) 
    102 #    plt.show() 
     117    #    plt.show() 
    103118     
    104 #    plt.figure(figsize=(12,5)) 
     119    #    plt.figure(figsize=(12,5)) 
    105120    c=ax2.contourf(xx,zz,ww,20) 
    106121    ax2.set_xlim((-.5,.5)), ax2.set_xlabel('x (km)') 
    … …  
    108123    plt.colorbar(c,ax=ax2) 
    109124    ax2.set_title('Vertical velocity at t=%g s (m/s)' % (it*T,)) 
    110 #    plt.tight_layout() 
    111 #    plt.show() 
     125    #    plt.tight_layout() 
     126    #    plt.show() 
    112127    plt.savefig('fig_NH_3D_bubble/%02d.png'%it) 
    113128     
    114129    time1, elapsed1 =time.time(), unst.getvar('elapsed') 
    115     flow = scheme.advance(flow,nt) 
     130    m,S,u,Phi,W = next_flow(m,S,u,Phi,W) 
    116131    time2, elapsed2 =time.time(), unst.getvar('elapsed') 
    117132    factor = 1000./(4*nt) 

  • codes/icosagcm/devel/Python/test/py/RSW2_MPAS_W02.py

    r641 r642  
    3838#scheme = time_step.RKn_simple(1,caldyn.bwd_fast_slow, dt) # forward Euler scheme 
    3939print dt, scheme.csjl, scheme.cfjl 
    40 step = unst.Caldyn_step(mesh,scheme) 
    41 step.setup_TRSW() 
     40step = unst.caldyn_step_TRSW(mesh,scheme,nt) 
    4241 
    4342u0 = Omega*radius/12. # cf Williamson (1991), p.13 
    … …  
    6362for i in range(20): 
    6463    if True: 
    65         for j in range(nt): 
    66             step.next() 
     64        step.next() 
    6765        gh,u = step.mass, step.u 
    6866        print i, gh.shape, gh.min(), gh.max(), u.min(), u.max() 

  • codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py

    r631 r642  
    1515    caldyn_thermo = unst.thermo_entropy 
    1616    g = mesh.dx/metric.dx_g0 
    17     caldyn = unst.Caldyn_HPE(caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g) 
    1817 
    1918    Sik = m0ik*gas0.s 
    … …  
    2524        S0ik = m0ik*gas0.s 
    2625 
    27     u0=mesh.field_u()  
     26    u0=mesh.field_u() 
    2827    u0[:,range(0,2*nx,2),:] = 100.*mesh.dx # Doppler shift 100 m/s 
    2928    flow0=(m0ik,S0ik,u0) 
    … …  
    3534    dt=T/nt 
    3635    print 'nt,dt,dx', nt,dt,dx 
    37     scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7) 
    38     #scheme = time_step.ARK3(caldyn.bwd_fast_slow, dt) 
    3936 
    40     flow=flow0 
     37    m,S,u=flow0 
     38 
     39    if False: # time stepping in Python 
     40        caldyn = unst.Caldyn_HPE(caldyn_thermo,caldyn_eta, mesh,thermo,BC,g) 
     41        scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7) 
     42        #scheme = time_step.ARK3(caldyn.bwd_fast_slow, dt) 
     43        def next_flow(m,S,u): 
     44            m,S,u = scheme.advance((m,S,u),nt) 
     45            return m,S,u,caldyn.geopot 
     46    else: # time stepping in Fortran 
     47        scheme = time_step.ARK2(None, dt, a32=0.7) 
     48        caldyn_step = unst.caldyn_step_HPE(mesh,scheme,nt, caldyn_thermo,caldyn_eta,thermo,BC,g) 
     49        def next_flow(m,S,u): 
     50            caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u 
     51            caldyn_step.next() 
     52            return caldyn_step.mass, caldyn_step.theta_rhodz, caldyn_step.u, caldyn_step.geopot 
     53         
    4154    for i in range(10): 
    4255        time1=time.time() 
    43         flow = scheme.advance(flow,nt) 
     56        m,S,u,geopot = next_flow(m,S,u) 
    4457        time2=time.time() 
    45         m,S,u = flow 
    46         print 'ms per time step : ', 1000*(time2-time1)/nt 
    47         print 'ptop, model top (m) :', unst.getvar('ptop'), caldyn.geopot.max()/unst.getvar('g') 
    48         junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) 
    49         zz=caldyn.geopot[:,0:llm]/metric.g0/1000 
     58        print 'ms per time step : ', 1000*(time2-time1)/nt, 1000*unst.getvar('elapsed')/(nt*(i+1)) 
     59        print 'ptop, model top (m) :', unst.getvar('ptop'), geopot.max()/unst.getvar('g') 
     60        #        junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) 
     61        zz=geopot[:,0:llm]/metric.g0/1000 
    5062        plt.figure(figsize=(12,3)) 
    5163        ucomp = mesh.ucomp(u) 

  • codes/icosagcm/devel/src/icosa_init.f90

    r583 r642  
    107107  END SUBROUTINE icosa_init 
    108108 
     109  SUBROUTINE force_recompile_unstructured 
     110    USE timestep_unstructured_mod 
     111  END SUBROUTINE force_recompile_unstructured 
     112 
    109113END MODULE icosa_init_mod 

  • codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90

    r638 r642  
    103103  FIELD_THETA :: dtheta_rhodz, theta 
    104104  FIELD_GEOPOT :: wflux 
    105   DOUBLE2(llm+1, edge_num) :: wwuu 
     105  FIELD_UL    :: wwuu 
    106106  DECLARE_INDICES 
    107107  DBL :: dF_deta, dFu_deta 

  • codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90

    r638 r642  
    3636#define HASNAN(field) (ANY(.NOT.ABS(field)<1e20)) 
    3737 
    38 SUBROUTINE caldyn_unstructured(tau, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! IN : flow state 
    39                                theta,ps,pk,hflux,qv, &    ! OUT : diags (except surface geopot : IN) 
    40                                dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, & 
    41                                dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(caldyn_unstructured) ! OUT : tendencies 
    42   DBL, VALUE :: tau 
    43   FIELD_MASS   :: rhodz, drhodz, pk, berni         ! IN, OUT, DIAG 
    44   FIELD_THETA  :: theta_rhodz, dtheta_rhodz, theta ! IN, OUT, DIAG 
    45   FIELD_GEOPOT :: wflux, w, geopot, &              ! DIAG, INOUT 
    46        dPhi_fast, dPhi_slow, dW_fast, dW_slow      ! OUT 
    47   FIELD_U      :: u,du_fast,du_slow,hflux,qu       ! INOUT,OUT,OUT,DIAG 
    48   FIELD_Z      :: qv                               ! DIAG 
    49   FIELD_PS     :: ps,dmass_col,mass_col            ! OUT,OUT,IN (if eta_mass) or OUT,UNUSED,UNUSED (if eta_lag) 
    50   DOUBLE2(llm+1, edge_num) :: wwuu 
    51   DBL          :: time1,time2 
    52   INTEGER :: ij 
    53  
    54 !  CALL CPU_TIME(time1) 
    55   time1=OMP_GET_WTIME() 
    56  
    57   IF(hydrostatic) THEN 
    58  
    59     !$OMP PARALLEL NUM_THREADS(nb_threads) 
    60     !$OMP DO SCHEDULE(STATIC) 
    61     DO ij=1,edge_num 
    62       du_fast(:,ij)=0. 
    63       du_slow(:,ij)=0. 
    64     END DO 
    65     !$OMP END DO 
    66     CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) 
    67     CALL compute_geopot(rhodz,theta, ps,pk,geopot) 
    68  
    69     CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u) 
    70  
    71     CALL compute_pvort_only(rhodz,u,qv,qu) 
    72     CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow) 
    73     CALL compute_coriolis(hflux,theta,qu, drhodz,dtheta_rhodz,du_slow) 
    74     IF(caldyn_eta == eta_mass) THEN 
    75        CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,wwuu) 
     38  SUBROUTINE caldyn_unstructured(tau, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! IN : flow state 
     39                                 theta,ps,pk,hflux,qv, &    ! OUT : diags (except surface geopot : IN) 
     40                                 dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, & 
     41                                 dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(caldyn_unstructured) ! OUT : tendencies 
     42    DBL, VALUE :: tau 
     43    FIELD_MASS   :: rhodz, drhodz, pk, berni         ! IN, OUT, DIAG 
     44    FIELD_THETA  :: theta_rhodz, dtheta_rhodz, theta ! IN, OUT, DIAG 
     45    FIELD_GEOPOT :: wflux, w, geopot, &              ! DIAG, INOUT 
     46         dPhi_fast, dPhi_slow, dW_fast, dW_slow      ! OUT 
     47    FIELD_U      :: u,du_fast,du_slow,hflux,qu       ! INOUT,OUT,OUT,DIAG 
     48    FIELD_Z      :: qv                               ! DIAG 
     49    FIELD_PS     :: ps,dmass_col,mass_col            ! OUT,OUT,IN (if eta_mass) or OUT,UNUSED,UNUSED (if eta_lag) 
     50    FIELD_UL     :: wwuu 
     51    DBL          :: time1,time2 
     52    INTEGER :: ij 
     53     
     54    !  CALL CPU_TIME(time1) 
     55    time1=OMP_GET_WTIME() 
     56     
     57    IF(hydrostatic) THEN 
     58        
     59       !$OMP PARALLEL NUM_THREADS(nb_threads) 
     60       !$OMP DO SCHEDULE(STATIC) 
     61       DO ij=1,edge_num 
     62          du_fast(:,ij)=0. 
     63          du_slow(:,ij)=0. 
     64       END DO 
     65       !$OMP END DO 
     66       CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) 
     67       CALL compute_geopot(rhodz,theta, ps,pk,geopot) 
     68        
     69       CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u) 
     70        
     71       CALL compute_pvort_only(rhodz,u,qv,qu) 
     72       CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow) 
     73       CALL compute_coriolis(hflux,theta,qu, drhodz,dtheta_rhodz,du_slow) 
     74       IF(caldyn_eta == eta_mass) THEN 
     75          CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,wwuu) 
     76       END IF 
     77       !$OMP END PARALLEL 
     78        
     79    ELSE ! NH 
     80       DO ij=1,edge_num 
     81          du_fast(:,ij)=0. 
     82          du_slow(:,ij)=0. 
     83       END DO 
     84       DO ij=1,primal_num 
     85          wflux(1,ij)=0. 
     86          wflux(llm+1,ij)=0. 
     87       END DO 
     88       CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) 
     89       CALL compute_caldyn_solver(tau,rhodz,theta,pk,geopot,W,dPhi_fast,dW_fast,du_fast) 
     90       CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u) 
     91       CALL compute_pvort_only(rhodz,u,qv,qu) 
     92       CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, hflux,du_slow,dPhi_slow,dW_slow)  
     93       CALL compute_coriolis(hflux,theta,qu, drhodz,dtheta_rhodz,du_slow) 
     94       IF(caldyn_eta == eta_mass) THEN 
     95          CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,wwuu) 
     96          CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, du_slow,dPhi_slow,dW_slow) 
     97       END IF 
    7698    END IF 
    77 !$OMP END PARALLEL 
    78  
    79  
    80   ELSE ! NH 
    81     DO ij=1,edge_num 
    82       du_fast(:,ij)=0. 
    83       du_slow(:,ij)=0. 
    84     END DO 
    85     DO ij=1,primal_num 
    86       wflux(1,ij)=0. 
    87       wflux(llm+1,ij)=0. 
    88     END DO 
    89     CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) 
    90     CALL compute_caldyn_solver(tau,rhodz,theta,pk,geopot,W,dPhi_fast,dW_fast,du_fast) 
    91     CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u) 
    92     CALL compute_pvort_only(rhodz,u,qv,qu) 
    93     CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, hflux,du_slow,dPhi_slow,dW_slow)  
    94     CALL compute_coriolis(hflux,theta,qu, drhodz,dtheta_rhodz,du_slow) 
    95     IF(caldyn_eta == eta_mass) THEN 
    96        CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,wwuu) 
    97        CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, du_slow,dPhi_slow,dW_slow) 
    98     END IF 
    99   END IF 
    100   
    101   time2=OMP_GET_WTIME() 
    102 !  CALL CPU_TIME(time2) 
    103   IF(time2>time1) elapsed = elapsed + time2-time1 
    104 END SUBROUTINE caldyn_unstructured 
    105 ! 
    106 !----------------------------- Time stepping ------------------------------- 
    107  
    108 ! 
    109   SUBROUTINE ARK_step(mass_col,rhodz,theta_rhodz,u,geopot,w, & ! INOUT : flow state 
     99     
     100    time2=OMP_GET_WTIME() 
     101    !  CALL CPU_TIME(time2) 
     102    IF(time2>time1) elapsed = elapsed + time2-time1 
     103  END SUBROUTINE caldyn_unstructured 
     104  ! 
     105  !----------------------------- Time stepping ------------------------------- 
     106   
     107  ! 
     108  SUBROUTINE ARK_step(nstep, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! INOUT : flow state 
    110109       theta,ps,pk,hflux,qv, &    ! OUT : diags (except surface geopot : IN) 
    111110       dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, & 
    112111       dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(ARK_step) ! OUT : tendencies 
     112    INTEGER, VALUE :: nstep ! advance by nstep time steps 
    113113    FIELD_PS     :: mass_col, ps                     ! OUT,IN (if eta_mass) or OUT,UNUSED (if eta_lag) 
    114114    FIELD_MASS   :: rhodz, pk, berni                 ! IN, DIAG 
    … …  
    123123    DOUBLE3(llm+1, primal_num, max_nb_stage) :: & 
    124124         dPhi_fast, dPhi_slow, dW_fast, dW_slow      ! OUT 
    125     FIELD_UL     :: DePhil, v_el, G_el               ! DIAG*3 
    126     DOUBLE2(llm+1, edge_num) :: wwuu 
     125    FIELD_UL     :: DePhil, v_el, G_el, wwuu         ! DIAG*4 
    127126    DBL       :: time1,time2 
    128     INTEGER :: stage, ij 
    129  
     127    INTEGER :: step, stage, ij 
     128     
    130129    !CALL CPU_TIME(time1) 
    131130    time1=OMP_GET_WTIME() 
    132  
     131     
    133132    !$OMP PARALLEL NUM_THREADS(nb_threads) 
    134133     
    135     DO stage=1, nb_stage 
    136         
    137        IF(hydrostatic) THEN 
    138  
    139           !$OMP DO SCHEDULE(STATIC) 
    140           DO ij=1,edge_num 
    141              du_fast(:,ij,stage)=0. 
    142              du_slow(:,ij,stage)=0. 
    143           END DO 
    144  
    145           CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) 
    146           CALL compute_geopot(rhodz,theta, ps,pk,geopot) 
    147            
    148           CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u) 
    149            
    150           CALL compute_pvort_only(rhodz,u,qv,qu) 
    151           CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage)) 
    152           CALL compute_coriolis(hflux,theta,qu, drhodz(:,:,stage), dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage)) 
    153           IF(caldyn_eta == eta_mass) THEN 
    154              STOP ! FIXME 
    155              CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, & 
    156                   dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),wwuu) 
     134    DO step=1, nstep 
     135       DO stage=1, nb_stage 
     136           
     137          IF(hydrostatic) THEN 
     138              
     139             !$OMP DO SCHEDULE(STATIC) 
     140             DO ij=1,edge_num 
     141                du_fast(:,ij,stage)=0. 
     142                du_slow(:,ij,stage)=0. 
     143             END DO 
     144              
     145             CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) 
     146             CALL compute_geopot(rhodz,theta, ps,pk,geopot) 
     147              
     148             CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u) 
     149              
     150             CALL compute_pvort_only(rhodz,u,qv,qu) 
     151             CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage)) 
     152             CALL compute_coriolis(hflux,theta,qu, drhodz(:,:,stage), dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage)) 
     153             IF(caldyn_eta == eta_mass) THEN 
     154                CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, & 
     155                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),wwuu) 
     156             END IF 
     157              
     158          ELSE ! NH 
     159             !$OMP DO SCHEDULE(STATIC) 
     160             DO ij=1,primal_num 
     161                wflux(1,ij)=0. 
     162                wflux(llm+1,ij)=0. 
     163             END DO 
     164             !$OMP END DO 
     165             CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) 
     166             CALL compute_caldyn_solver(tauj(stage),rhodz,theta,pk,geopot,W, & 
     167                  dPhi_fast(:,:,stage), dW_fast(:,:,stage), du_fast(:,:,stage)) 
     168             CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage),u) 
     169             CALL compute_pvort_only(rhodz,u,qv,qu) 
     170             CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, hflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage)) 
     171             CALL compute_coriolis(hflux,theta,qu, drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage)) 
     172             IF(caldyn_eta == eta_mass) THEN 
     173                CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, & 
     174                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), wwuu) 
     175                CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) ) 
     176             END IF 
     177          END IF ! NH 
     178           
     179          ! FIXME : mass_col is computed from rhodz 
     180          ! so that the DOFs are the same whatever caldyn_eta 
     181          ! in DYNAMICO mass_col is prognosed rather than rhodz 
     182           
     183#define UPDATE(clj,field,dfield) update(stage,SIZE(field),clj,field,dfield) 
     184           
     185          CALL UPDATE(cslj, rhodz, drhodz) 
     186          CALL UPDATE(cslj, theta_rhodz, dtheta_rhodz) 
     187          CALL UPDATE(cslj, u, du_slow) 
     188          CALL UPDATE(cflj, u, du_fast) 
     189           
     190          IF(.NOT.hydrostatic) THEN 
     191             CALL UPDATE(cslj, geopot, dPhi_slow) 
     192             CALL UPDATE(cflj, geopot, dPhi_fast) 
     193             CALL UPDATE(cslj, W, dW_slow) 
     194             CALL UPDATE(cflj, W, dW_fast) 
    157195          END IF 
    158196           
    159        ELSE ! NH 
    160           STOP ! FIXME 
    161           !$OMP DO SCHEDULE(STATIC) 
    162           DO ij=1,primal_num 
    163              wflux(1,ij)=0. 
    164              wflux(llm+1,ij)=0. 
    165           END DO 
    166           !$OMP END DO 
    167           CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) 
    168           CALL compute_caldyn_solver(tauj(stage),rhodz,theta,pk,geopot,W, & 
    169                dPhi_fast(:,:,stage), dW_fast(:,:,stage), du_fast(:,:,stage)) 
    170           CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage),u) 
    171           CALL compute_pvort_only(rhodz,u,qv,qu) 
    172           CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, hflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage)) 
    173           CALL compute_coriolis(hflux,theta,qu, drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage)) 
    174           IF(caldyn_eta == eta_mass) THEN 
    175              CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, & 
    176                   dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), wwuu) 
    177              CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) ) 
    178           END IF 
    179        END IF ! NH 
    180         
    181       ! FIXME : mass_col is computed from rhodz 
    182       ! so that the DOFs are the same whatever caldyn_eta 
    183       ! in DYNAMICO mass_col is prognosed rather than rhodz 
    184  
    185 #define UPDATE(clj,field,dfield) update(stage,SIZE(field),clj,field,dfield) 
    186  
    187        CALL UPDATE(cslj, rhodz, drhodz) 
    188        CALL UPDATE(cslj, theta_rhodz, dtheta_rhodz) 
    189        CALL UPDATE(cslj, u, du_slow) 
    190        CALL UPDATE(cflj, u, du_fast) 
    191  
    192        IF(.NOT.hydrostatic) THEN 
    193           STOP ! FIXME 
    194           CALL UPDATE(cslj, geopot, dPhi_slow) 
    195           CALL UPDATE(cflj, geopot, dPhi_fast) 
    196           CALL UPDATE(cslj, W, dW_slow) 
    197           CALL UPDATE(cflj, W, dW_fast) 
    198        END IF 
    199  
     197       END DO 
    200198    END DO 
    201 !$OMP END PARALLEL 
    202      
    203   time2=OMP_GET_WTIME() 
    204 !  CALL CPU_TIME(time2) 
    205   IF(time2>time1) elapsed = elapsed + time2-time1 
    206  
     199    !$OMP END PARALLEL 
     200     
     201    time2=OMP_GET_WTIME() 
     202    !  CALL CPU_TIME(time2) 
     203    IF(time2>time1) elapsed = elapsed + time2-time1 
     204     
    207205  END SUBROUTINE ARK_step 
    208   ! 
    209 ! 
    210 SUBROUTINE update(j,sz,clj,field,dfield) 
    211   INTEGER :: j, sz ! stage in ARK scheme, field size 
    212   DOUBLE2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau 
    213   DOUBLE1(sz) :: field 
    214   DOUBLE2(sz, max_nb_stage) :: dfield 
    215   ! 
    216   INTEGER :: l, ij 
    217 !  PRINT *, clj(1:j,j) 
    218   SELECT CASE(j) 
    219   CASE(1) 
    220      !$OMP DO SCHEDULE(static) 
    221      DO ij=1,sz 
    222         field(ij) = field(ij) & 
    223              + clj(1,j)*dfield(ij,1) 
    224      END DO 
    225   CASE(2) 
    226      !$OMP DO SCHEDULE(static) 
    227      DO ij=1,sz 
    228         field(ij) = field(ij) & 
    229              + clj(1,j)*dfield(ij,1) & 
    230              + clj(2,j)*dfield(ij,2) 
    231      END DO 
    232   CASE(3) 
    233      !$OMP DO SCHEDULE(static) 
    234      DO ij=1,sz 
    235         field(ij) = field(ij) & 
    236              + clj(1,j)*dfield(ij,1) & 
    237              + clj(2,j)*dfield(ij,2) & 
    238              + clj(3,j)*dfield(ij,3) 
    239  
    240      END DO 
    241   CASE(4) 
    242      !$OMP DO SCHEDULE(static) 
    243      DO ij=1,sz 
    244         field(ij) = field(ij) & 
    245              + clj(1,j)*dfield(ij,1) & 
    246              + clj(2,j)*dfield(ij,2) & 
    247              + clj(3,j)*dfield(ij,3) & 
    248              + clj(4,j)*dfield(ij,4) 
    249      END DO 
    250   CASE(5) 
    251      !$OMP DO SCHEDULE(static) 
    252      DO ij=1,sz 
    253         field(ij) = field(ij) & 
    254              + clj(1,j)*dfield(ij,1) & 
    255              + clj(2,j)*dfield(ij,2) & 
    256              + clj(3,j)*dfield(ij,3) & 
    257              + clj(4,j)*dfield(ij,4) & 
    258              + clj(5,j)*dfield(ij,5) 
    259      END DO 
    260   END SELECT 
    261 END SUBROUTINE update 
    262  
    263 !---------------------------------------------- XIOS ---------------------------------------- 
     206   
     207   
     208  SUBROUTINE update(j,sz,clj,field,dfield) 
     209    INTEGER :: j, sz ! stage in ARK scheme, field size 
     210    DOUBLE2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau 
     211    DOUBLE1(sz) :: field 
     212    DOUBLE2(sz, max_nb_stage) :: dfield 
     213    ! 
     214    INTEGER :: l, ij 
     215    !  PRINT *, clj(1:j,j) 
     216    SELECT CASE(j) 
     217    CASE(1) 
     218       !$OMP DO SCHEDULE(static) 
     219       DO ij=1,sz 
     220          field(ij) = field(ij) & 
     221               + clj(1,j)*dfield(ij,1) 
     222       END DO 
     223    CASE(2) 
     224       !$OMP DO SCHEDULE(static) 
     225       DO ij=1,sz 
     226          field(ij) = field(ij) & 
     227               + clj(1,j)*dfield(ij,1) & 
     228               + clj(2,j)*dfield(ij,2) 
     229       END DO 
     230    CASE(3) 
     231       !$OMP DO SCHEDULE(static) 
     232       DO ij=1,sz 
     233          field(ij) = field(ij) & 
     234               + clj(1,j)*dfield(ij,1) & 
     235               + clj(2,j)*dfield(ij,2) & 
     236               + clj(3,j)*dfield(ij,3) 
     237           
     238       END DO 
     239    CASE(4) 
     240       !$OMP DO SCHEDULE(static) 
     241       DO ij=1,sz 
     242          field(ij) = field(ij) & 
     243               + clj(1,j)*dfield(ij,1) & 
     244               + clj(2,j)*dfield(ij,2) & 
     245               + clj(3,j)*dfield(ij,3) & 
     246               + clj(4,j)*dfield(ij,4) 
     247       END DO 
     248    CASE(5) 
     249       !$OMP DO SCHEDULE(static) 
     250       DO ij=1,sz 
     251          field(ij) = field(ij) & 
     252               + clj(1,j)*dfield(ij,1) & 
     253               + clj(2,j)*dfield(ij,2) & 
     254               + clj(3,j)*dfield(ij,3) & 
     255               + clj(4,j)*dfield(ij,4) & 
     256               + clj(5,j)*dfield(ij,5) 
     257       END DO 
     258    END SELECT 
     259  END SUBROUTINE update 
     260   
     261  !---------------------------------------------- XIOS ---------------------------------------- 
    264262 
    265263#ifdef CPP_USING_XIOS 
    266  
     264   
    267265  SUBROUTINE setup_xios() BINDC(setup_xios) 
    268266    ! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine 
Note: See TracChangeset for help on using the changeset viewer.