Changeset 642 for codes/icosagcm/devel
- Timestamp:
- 12/19/17 15:26:51 (7 years ago)
- Location:
- codes/icosagcm/devel
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/dynamico/meshes.py
r639 r642 148 148 Aiv,Aiv,f+0.*Aiv,le_de,Riv2,-wee) 149 149 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)) 156 160 def ucomp(self,u): return u[:,range(0,2*self.nx,2),:] 157 161 def set_ucomp(self,uv,u): uv[:,range(0,2*self.nx,2),:]=u -
codes/icosagcm/devel/Python/src/functions.h
r639 r642 7 7 void dynamico_init_params(void); 8 8 9 void dynamico_ARK_step(double *mass_col, double *rhodz, double *theta_rhodz, 9 void dynamico_ARK_step(int nstep, 10 double *mass_col, double *rhodz, double *theta_rhodz, 10 11 double *u, double *geopot, double *w, 11 12 double *theta, double *ps, double *pk, double *hflux, double *qv, -
codes/icosagcm/devel/Python/src/unstructured.pyx
r640 r642 18 18 19 19 cdef 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, 21 22 double *u, double *geopot, double *w, 22 23 double *theta, double *ps, double *pk, double *hflux, double *qv, … … 78 79 cdef dbl_ptr ptr2(double[:,:] data) : return &data[0,0] 79 80 cdef dbl_ptr ptr3(double[:,:,:] data) : return &data[0,0,0] 81 cdef dbl_ptr ptr4(double[:,:,:,:] data) : return &data[0,0,0,0] 80 82 cdef dbl_ptr ptr(data): 81 83 n=data.ndim … … 83 85 if n==2 : return ptr2(data) 84 86 if n==3 : return ptr3(data) 85 87 if n==4 : return ptr4(data) 88 if n>4: raise IndexError 89 86 90 cdef alloc(dbl_ptr *p, allocator, n=1): 87 91 data=allocator(n) … … 93 97 94 98 cdef class Caldyn_step: 99 # number of time steps to do at each invocation of advance() 100 cdef int nstep 95 101 # pointer to allocated arrays 96 cdef dbl_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_ w# prognostic102 cdef dbl_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_W # prognostic 97 103 cdef dbl_ptr p_mass_col, p_dmass_col, p_ps, p_theta, p_pk, p_hflux, p_qv # diagnostic 98 104 cdef dbl_ptr p_drhodz, p_dtheta_rhodz, p_du_fast, p_du_slow # tendencies 99 105 cdef dbl_ptr p_dPhi_fast, p_dPhi_slow, p_dW_fast, p_dW_slow # tendencies 100 106 # allocated arrays, must remain referenced or segfault 101 cdef readonly np.ndarray mass, theta_rhodz, u, geopot, w107 cdef readonly np.ndarray mass, theta_rhodz, u, geopot, W 102 108 cdef readonly np.ndarray mass_col, dmass_col, ps, theta, pk, hflux, qv 103 109 cdef readonly np.ndarray drhodz, dtheta_rhodz, du_fast, du_slow 104 110 cdef readonly np.ndarray dPhi_fast, dPhi_slow, dW_fast, dW_slow 105 111 106 def __init__(self,mesh,time_scheme): 112 def __init__(self,mesh,time_scheme, nstep): 113 self.nstep=nstep 107 114 # self.mesh=mesh 108 115 fps, ftheta, fmass = mesh.field_ps, mesh.field_theta, mesh.field_mass … … 113 120 cdef double[:,:] cflj_ = time_scheme.cfjl 114 121 ns = time_scheme.nstage 115 nb_stage[0]= ns 122 nb_stage[0]=ns 123 116 124 cdef int i,j 117 125 for i in range(ns): … … 123 131 # prognostic/diagnostic 124 132 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), 126 134 self.mass, self.theta_rhodz = alloc(&self.p_mass, fmass), alloc(&self.p_theta_rhodz, fmass), 127 135 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), 129 137 self.hflux, self.u = alloc(&self.p_hflux, fu), alloc(&self.p_u, fu) 130 138 self.qv = alloc(&self.p_qv,fz) … … 137 145 # global elapsed 138 146 # 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, 143 150 self.p_theta, self.p_ps, self.p_pk, self.p_hflux, self.p_qv, 144 151 self.p_dmass_col, self.p_drhodz, self.p_dtheta_rhodz, … … 148 155 #time2=time.time() 149 156 #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 158 def 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) 163 def 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) 170 def 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 169 179 #----------------------------- Base class for dynamics ------------------------ 170 180 … … 212 222 213 223 class 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): 215 225 Caldyn.__init__(self,mesh) 216 226 setvars(('hydrostatic','caldyn_thermo','caldyn_eta', … … 231 241 232 242 class 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): 234 244 Caldyn.__init__(self,mesh) 235 245 setvars(('hydrostatic','caldyn_thermo','caldyn_eta', -
codes/icosagcm/devel/Python/test/fig_RSW2_MPAS_W02
-
Property
svn:ignore
set to
*.png
-
Property
svn:ignore
set to
-
codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py
r632 r642 53 53 return thermo, mesh, params, (mass_ik,Sik,ujk,Phi_il,Wil), gas 54 54 55 Lx, nx, llm, thetac, T, Nslice, courant = 2000., 100, 50, 30., 5., 10, 2.856 #Lx, nx, llm, thetac, T, Nslice, courant = 2000., 50, 25, 1., 50., 10, 2.855 #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, 30, 5., 10, 2.8 57 57 #Lx, nx, llm, thetac, T, Nslice, courant = 3000., 75, 25, -30, 5., 10, 2.8 58 58 … … 61 61 thermo, mesh, params, flow0, gas0 = thermal_bubble_3D(Lx,nx,Ly,ny,llm, thetac=thetac) 62 62 63 caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass64 caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,params,thermo,params,params.g)65 66 63 # compute hybrid coefs from initial distribution of mass 67 64 mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0]) 68 65 unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk) 69 print mass_dbk70 66 71 67 dz = flow0[3].max()/(params.g*llm) … … 74 70 nt = int(math.ceil(T/dt)) 75 71 dt = T/nt 76 print 'Time step : % g s' % dt72 print 'Time step : %d x %g s' % (nt,dt) 77 73 78 scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt) 74 caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass 75 #caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_lag 79 76 80 flow=flow0 77 if 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 84 else: # 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 96 m,S,u,Phi,W=flow0 81 97 w=mesh.field_mass() 82 98 z=mesh.field_mass() 99 83 100 for 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)) 87 102 for l in range(llm): 88 103 w[:,:,l]=.5*params.g*(W[:,:,l+1]+W[:,:,l])/m[:,:,l] … … 92 107 jj=ny/2 93 108 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)) 94 111 95 f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12,4))96 97 112 c=ax1.contourf(xx,zz,ss,20) 98 113 ax1.set_xlim((-.5,.5)), ax1.set_xlabel('x (km)') … … 100 115 plt.colorbar(c,ax=ax1) 101 116 ax1.set_title('Specific entropy at t=%g s (J/K/kg)' % (it*T,)) 102 # plt.show()117 # plt.show() 103 118 104 # plt.figure(figsize=(12,5))119 # plt.figure(figsize=(12,5)) 105 120 c=ax2.contourf(xx,zz,ww,20) 106 121 ax2.set_xlim((-.5,.5)), ax2.set_xlabel('x (km)') … … 108 123 plt.colorbar(c,ax=ax2) 109 124 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() 112 127 plt.savefig('fig_NH_3D_bubble/%02d.png'%it) 113 128 114 129 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) 116 131 time2, elapsed2 =time.time(), unst.getvar('elapsed') 117 132 factor = 1000./(4*nt) -
codes/icosagcm/devel/Python/test/py/RSW2_MPAS_W02.py
r641 r642 38 38 #scheme = time_step.RKn_simple(1,caldyn.bwd_fast_slow, dt) # forward Euler scheme 39 39 print dt, scheme.csjl, scheme.cfjl 40 step = unst.Caldyn_step(mesh,scheme) 41 step.setup_TRSW() 40 step = unst.caldyn_step_TRSW(mesh,scheme,nt) 42 41 43 42 u0 = Omega*radius/12. # cf Williamson (1991), p.13 … … 63 62 for i in range(20): 64 63 if True: 65 for j in range(nt): 66 step.next() 64 step.next() 67 65 gh,u = step.mass, step.u 68 66 print i, gh.shape, gh.min(), gh.max(), u.min(), u.max() -
codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py
r631 r642 15 15 caldyn_thermo = unst.thermo_entropy 16 16 g = mesh.dx/metric.dx_g0 17 caldyn = unst.Caldyn_HPE(caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g)18 17 19 18 Sik = m0ik*gas0.s … … 25 24 S0ik = m0ik*gas0.s 26 25 27 u0=mesh.field_u() 26 u0=mesh.field_u() 28 27 u0[:,range(0,2*nx,2),:] = 100.*mesh.dx # Doppler shift 100 m/s 29 28 flow0=(m0ik,S0ik,u0) … … 35 34 dt=T/nt 36 35 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)39 36 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 41 54 for i in range(10): 42 55 time1=time.time() 43 flow = scheme.advance(flow,nt)56 m,S,u,geopot = next_flow(m,S,u) 44 57 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 50 62 plt.figure(figsize=(12,3)) 51 63 ucomp = mesh.ucomp(u) -
codes/icosagcm/devel/src/icosa_init.f90
r583 r642 107 107 END SUBROUTINE icosa_init 108 108 109 SUBROUTINE force_recompile_unstructured 110 USE timestep_unstructured_mod 111 END SUBROUTINE force_recompile_unstructured 112 109 113 END MODULE icosa_init_mod -
codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90
r638 r642 103 103 FIELD_THETA :: dtheta_rhodz, theta 104 104 FIELD_GEOPOT :: wflux 105 DOUBLE2(llm+1, edge_num):: wwuu105 FIELD_UL :: wwuu 106 106 DECLARE_INDICES 107 107 DBL :: dF_deta, dFu_deta -
codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90
r638 r642 36 36 #define HASNAN(field) (ANY(.NOT.ABS(field)<1e20)) 37 37 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 76 98 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 110 109 theta,ps,pk,hflux,qv, & ! OUT : diags (except surface geopot : IN) 111 110 dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, & 112 111 dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(ARK_step) ! OUT : tendencies 112 INTEGER, VALUE :: nstep ! advance by nstep time steps 113 113 FIELD_PS :: mass_col, ps ! OUT,IN (if eta_mass) or OUT,UNUSED (if eta_lag) 114 114 FIELD_MASS :: rhodz, pk, berni ! IN, DIAG … … 123 123 DOUBLE3(llm+1, primal_num, max_nb_stage) :: & 124 124 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 127 126 DBL :: time1,time2 128 INTEGER :: st age, ij129 127 INTEGER :: step, stage, ij 128 130 129 !CALL CPU_TIME(time1) 131 130 time1=OMP_GET_WTIME() 132 131 133 132 !$OMP PARALLEL NUM_THREADS(nb_threads) 134 133 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) 157 195 END IF 158 196 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 200 198 END DO 201 !$OMP END PARALLEL202 203 time2=OMP_GET_WTIME()204 ! CALL CPU_TIME(time2)205 IF(time2>time1) elapsed = elapsed + time2-time1206 199 !$OMP END PARALLEL 200 201 time2=OMP_GET_WTIME() 202 ! CALL CPU_TIME(time2) 203 IF(time2>time1) elapsed = elapsed + time2-time1 204 207 205 END SUBROUTINE ARK_step 208 !209 ! 210 SUBROUTINE update(j,sz,clj,field,dfield)211 INTEGER :: j, sz ! stage in ARK scheme, field size212 DOUBLE2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau213 DOUBLE1(sz) :: field214 DOUBLE2(sz, max_nb_stage) :: dfield215 !216 INTEGER :: l, ij217 ! PRINT *, clj(1:j,j)218 SELECT CASE(j)219 CASE(1)220 !$OMP DO SCHEDULE(static)221 DO ij=1,sz222 field(ij) = field(ij) &223 + clj(1,j)*dfield(ij,1)224 END DO225 CASE(2)226 !$OMP DO SCHEDULE(static)227 DO ij=1,sz228 field(ij) = field(ij) &229 + clj(1,j)*dfield(ij,1) &230 + clj(2,j)*dfield(ij,2)231 END DO232 CASE(3)233 !$OMP DO SCHEDULE(static)234 DO ij=1,sz235 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 DO241 CASE(4)242 !$OMP DO SCHEDULE(static)243 DO ij=1,sz244 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 DO250 CASE(5)251 !$OMP DO SCHEDULE(static)252 DO ij=1,sz253 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 DO260 END SELECT261 END SUBROUTINE update262 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 ---------------------------------------- 264 262 265 263 #ifdef CPP_USING_XIOS 266 264 267 265 SUBROUTINE setup_xios() BINDC(setup_xios) 268 266 ! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine
Note: See TracChangeset
for help on using the changeset viewer.