[638] | 1 | MODULE timestep_unstructured_mod |
---|
| 2 | USE ISO_C_BINDING |
---|
| 3 | USE OMP_LIB |
---|
| 4 | #ifdef CPP_USING_XIOS |
---|
| 5 | USE xios |
---|
| 6 | #endif |
---|
| 7 | USE caldyn_unstructured_mod |
---|
| 8 | IMPLICIT NONE |
---|
| 9 | PRIVATE |
---|
| 10 | SAVE |
---|
| 11 | |
---|
| 12 | #define BINDC_(thename) BIND(C, name=#thename) |
---|
| 13 | #define BINDC(thename) BINDC_(dynamico_ ## thename) |
---|
| 14 | |
---|
| 15 | #define DBL REAL(C_DOUBLE) |
---|
| 16 | #define DOUBLE1(m) DBL, DIMENSION(m) |
---|
| 17 | #define DOUBLE2(m,n) DBL, DIMENSION(m,n) |
---|
| 18 | #define DOUBLE3(m,n,p) DBL, DIMENSION(m,n,p) |
---|
| 19 | #define DOUBLE4(m,n,p,q) DBL, DIMENSION(m,n,p,q) |
---|
| 20 | #define INDEX INTEGER(C_INT) |
---|
| 21 | |
---|
| 22 | #ifdef CPP_USING_XIOS |
---|
| 23 | TYPE(xios_context) :: ctx_hdl |
---|
| 24 | #endif |
---|
| 25 | |
---|
| 26 | CONTAINS |
---|
| 27 | |
---|
| 28 | #define FIELD_PS DOUBLE1(primal_num) |
---|
| 29 | #define FIELD_MASS DOUBLE2(llm, primal_num) |
---|
| 30 | #define FIELD_Z DOUBLE2(llm, dual_num) |
---|
| 31 | #define FIELD_U DOUBLE2(llm, edge_num) |
---|
| 32 | #define FIELD_UL DOUBLE2(llm+1, edge_num) |
---|
| 33 | #define FIELD_THETA DOUBLE3(llm, primal_num, nqdyn) |
---|
| 34 | #define FIELD_GEOPOT DOUBLE2(llm+1, primal_num) |
---|
| 35 | |
---|
| 36 | #define HASNAN(field) (ANY(.NOT.ABS(field)<1e20)) |
---|
| 37 | |
---|
[642] | 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 |
---|
[645] | 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 |
---|
[665] | 47 | FIELD_U :: u,du_fast,du_slow,hflux,qu ! INOUT,OUT,OUT,DIAG |
---|
[645] | 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) |
---|
[665] | 50 | ! buffers for fields that need to be shared among OpenMP threads |
---|
| 51 | ! vertical size is allocated as llm+1 even if only llm is needed |
---|
| 52 | FIELD_GEOPOT :: bufm1, bufm2, bufm3 |
---|
| 53 | FIELD_UL :: bufu1, bufu2, bufu3,bufu4 |
---|
| 54 | |
---|
[642] | 55 | DBL :: time1,time2 |
---|
| 56 | INTEGER :: ij |
---|
| 57 | |
---|
| 58 | ! CALL CPU_TIME(time1) |
---|
| 59 | time1=OMP_GET_WTIME() |
---|
| 60 | |
---|
| 61 | IF(hydrostatic) THEN |
---|
| 62 | |
---|
| 63 | !$OMP PARALLEL NUM_THREADS(nb_threads) |
---|
| 64 | !$OMP DO SCHEDULE(STATIC) |
---|
| 65 | DO ij=1,edge_num |
---|
| 66 | du_fast(:,ij)=0. |
---|
| 67 | du_slow(:,ij)=0. |
---|
| 68 | END DO |
---|
| 69 | !$OMP END DO |
---|
| 70 | CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) |
---|
| 71 | CALL compute_geopot(rhodz,theta, ps,pk,geopot) |
---|
| 72 | |
---|
| 73 | CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u) |
---|
| 74 | |
---|
| 75 | CALL compute_pvort_only(rhodz,u,qv,qu) |
---|
| 76 | CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow) |
---|
[665] | 77 | CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow) |
---|
[642] | 78 | IF(caldyn_eta == eta_mass) THEN |
---|
[665] | 79 | CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1) |
---|
[642] | 80 | END IF |
---|
| 81 | !$OMP END PARALLEL |
---|
| 82 | |
---|
| 83 | ELSE ! NH |
---|
| 84 | DO ij=1,edge_num |
---|
| 85 | du_fast(:,ij)=0. |
---|
| 86 | du_slow(:,ij)=0. |
---|
| 87 | END DO |
---|
| 88 | DO ij=1,primal_num |
---|
| 89 | wflux(1,ij)=0. |
---|
| 90 | wflux(llm+1,ij)=0. |
---|
| 91 | END DO |
---|
| 92 | CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) |
---|
[665] | 93 | CALL compute_caldyn_solver(tau,rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W,dPhi_fast,dW_fast,du_fast) |
---|
[642] | 94 | CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u) |
---|
| 95 | CALL compute_pvort_only(rhodz,u,qv,qu) |
---|
[665] | 96 | CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3,bufu1,bufu2,bufu3,bufu4, hflux,du_slow,dPhi_slow,dW_slow) |
---|
| 97 | CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow) |
---|
[642] | 98 | IF(caldyn_eta == eta_mass) THEN |
---|
[665] | 99 | CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1) |
---|
| 100 | CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, du_slow,dPhi_slow,dW_slow) |
---|
[642] | 101 | END IF |
---|
[638] | 102 | END IF |
---|
[642] | 103 | |
---|
| 104 | time2=OMP_GET_WTIME() |
---|
| 105 | ! CALL CPU_TIME(time2) |
---|
| 106 | IF(time2>time1) elapsed = elapsed + time2-time1 |
---|
[652] | 107 | |
---|
| 108 | CALL print_trace() |
---|
[642] | 109 | END SUBROUTINE caldyn_unstructured |
---|
| 110 | ! |
---|
| 111 | !----------------------------- Time stepping ------------------------------- |
---|
| 112 | |
---|
| 113 | ! |
---|
| 114 | SUBROUTINE ARK_step(nstep, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! INOUT : flow state |
---|
[638] | 115 | theta,ps,pk,hflux,qv, & ! OUT : diags (except surface geopot : IN) |
---|
| 116 | dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, & |
---|
| 117 | dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(ARK_step) ! OUT : tendencies |
---|
[642] | 118 | INTEGER, VALUE :: nstep ! advance by nstep time steps |
---|
[638] | 119 | FIELD_PS :: mass_col, ps ! OUT,IN (if eta_mass) or OUT,UNUSED (if eta_lag) |
---|
| 120 | FIELD_MASS :: rhodz, pk, berni ! IN, DIAG |
---|
| 121 | FIELD_THETA :: theta_rhodz, theta ! IN, DIAG |
---|
[665] | 122 | FIELD_U :: u,hflux,qu ! INOUT,DIAG*2 |
---|
| 123 | FIELD_GEOPOT :: geopot, w, wflux ! IN, INOUT, DIAG |
---|
[638] | 124 | FIELD_Z :: qv ! DIAG |
---|
| 125 | DOUBLE2( primal_num, max_nb_stage) :: dmass_col ! OUT |
---|
| 126 | DOUBLE3(llm, primal_num, max_nb_stage) :: drhodz ! OUT |
---|
| 127 | DOUBLE4(llm, primal_num, nqdyn, max_nb_stage) :: dtheta_rhodz ! OUT |
---|
| 128 | DOUBLE3(llm, edge_num, max_nb_stage) :: du_fast,du_slow ! OUT |
---|
| 129 | DOUBLE3(llm+1, primal_num, max_nb_stage) :: & |
---|
| 130 | dPhi_fast, dPhi_slow, dW_fast, dW_slow ! OUT |
---|
[665] | 131 | ! buffers for fields that need to be shared among OpenMP threads |
---|
| 132 | ! vertical size is allocated as llm+1 even if only llm is needed |
---|
| 133 | FIELD_GEOPOT :: bufm1, bufm2, bufm3 |
---|
| 134 | FIELD_UL :: bufu1, bufu2, bufu3,bufu4 |
---|
[638] | 135 | DBL :: time1,time2 |
---|
[642] | 136 | INTEGER :: step, stage, ij |
---|
| 137 | |
---|
[638] | 138 | !CALL CPU_TIME(time1) |
---|
| 139 | time1=OMP_GET_WTIME() |
---|
[642] | 140 | |
---|
[638] | 141 | !$OMP PARALLEL NUM_THREADS(nb_threads) |
---|
| 142 | |
---|
[642] | 143 | DO step=1, nstep |
---|
| 144 | DO stage=1, nb_stage |
---|
[638] | 145 | |
---|
[642] | 146 | IF(hydrostatic) THEN |
---|
| 147 | |
---|
| 148 | !$OMP DO SCHEDULE(STATIC) |
---|
| 149 | DO ij=1,edge_num |
---|
| 150 | du_fast(:,ij,stage)=0. |
---|
| 151 | du_slow(:,ij,stage)=0. |
---|
| 152 | END DO |
---|
| 153 | |
---|
| 154 | CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) |
---|
| 155 | CALL compute_geopot(rhodz,theta, ps,pk,geopot) |
---|
| 156 | |
---|
| 157 | CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u) |
---|
| 158 | |
---|
| 159 | CALL compute_pvort_only(rhodz,u,qv,qu) |
---|
| 160 | CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage)) |
---|
[665] | 161 | CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz(:,:,stage), dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage)) |
---|
[642] | 162 | IF(caldyn_eta == eta_mass) THEN |
---|
| 163 | CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, & |
---|
[665] | 164 | dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),bufu1) |
---|
[642] | 165 | END IF |
---|
| 166 | |
---|
| 167 | ELSE ! NH |
---|
| 168 | !$OMP DO SCHEDULE(STATIC) |
---|
| 169 | DO ij=1,primal_num |
---|
| 170 | wflux(1,ij)=0. |
---|
| 171 | wflux(llm+1,ij)=0. |
---|
| 172 | END DO |
---|
| 173 | !$OMP END DO |
---|
| 174 | CALL compute_theta(mass_col,rhodz,theta_rhodz, theta) |
---|
[665] | 175 | CALL compute_caldyn_solver(tauj(stage),rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W, & |
---|
[642] | 176 | dPhi_fast(:,:,stage), dW_fast(:,:,stage), du_fast(:,:,stage)) |
---|
| 177 | CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage),u) |
---|
| 178 | CALL compute_pvort_only(rhodz,u,qv,qu) |
---|
[665] | 179 | CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3,bufu1,bufu2,bufu3,bufu4, hflux, & |
---|
| 180 | du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage)) |
---|
| 181 | CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage)) |
---|
[642] | 182 | IF(caldyn_eta == eta_mass) THEN |
---|
| 183 | CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, & |
---|
[665] | 184 | dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), bufu1) |
---|
| 185 | CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, & |
---|
| 186 | du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) ) |
---|
[642] | 187 | END IF |
---|
| 188 | END IF ! NH |
---|
[638] | 189 | |
---|
[642] | 190 | ! FIXME : mass_col is computed from rhodz |
---|
| 191 | ! so that the DOFs are the same whatever caldyn_eta |
---|
| 192 | ! in DYNAMICO mass_col is prognosed rather than rhodz |
---|
| 193 | |
---|
| 194 | #define UPDATE(clj,field,dfield) update(stage,SIZE(field),clj,field,dfield) |
---|
| 195 | |
---|
| 196 | CALL UPDATE(cslj, rhodz, drhodz) |
---|
| 197 | CALL UPDATE(cslj, theta_rhodz, dtheta_rhodz) |
---|
| 198 | CALL UPDATE(cslj, u, du_slow) |
---|
| 199 | CALL UPDATE(cflj, u, du_fast) |
---|
| 200 | |
---|
| 201 | IF(.NOT.hydrostatic) THEN |
---|
| 202 | CALL UPDATE(cslj, geopot, dPhi_slow) |
---|
| 203 | CALL UPDATE(cflj, geopot, dPhi_fast) |
---|
| 204 | CALL UPDATE(cslj, W, dW_slow) |
---|
| 205 | CALL UPDATE(cflj, W, dW_fast) |
---|
[638] | 206 | END IF |
---|
| 207 | |
---|
[642] | 208 | END DO |
---|
[638] | 209 | END DO |
---|
[642] | 210 | !$OMP END PARALLEL |
---|
[638] | 211 | |
---|
[642] | 212 | time2=OMP_GET_WTIME() |
---|
| 213 | ! CALL CPU_TIME(time2) |
---|
| 214 | IF(time2>time1) elapsed = elapsed + time2-time1 |
---|
[652] | 215 | |
---|
| 216 | CALL print_trace() |
---|
[642] | 217 | |
---|
[638] | 218 | END SUBROUTINE ARK_step |
---|
[642] | 219 | |
---|
| 220 | |
---|
| 221 | SUBROUTINE update(j,sz,clj,field,dfield) |
---|
| 222 | INTEGER :: j, sz ! stage in ARK scheme, field size |
---|
| 223 | DOUBLE2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau |
---|
| 224 | DOUBLE1(sz) :: field |
---|
| 225 | DOUBLE2(sz, max_nb_stage) :: dfield |
---|
| 226 | ! |
---|
| 227 | INTEGER :: l, ij |
---|
[658] | 228 | CALL enter_trace(id_update, 8*sz*(j+1) ) |
---|
[642] | 229 | ! PRINT *, clj(1:j,j) |
---|
| 230 | SELECT CASE(j) |
---|
| 231 | CASE(1) |
---|
| 232 | !$OMP DO SCHEDULE(static) |
---|
| 233 | DO ij=1,sz |
---|
| 234 | field(ij) = field(ij) & |
---|
| 235 | + clj(1,j)*dfield(ij,1) |
---|
| 236 | END DO |
---|
| 237 | CASE(2) |
---|
| 238 | !$OMP DO SCHEDULE(static) |
---|
| 239 | DO ij=1,sz |
---|
| 240 | field(ij) = field(ij) & |
---|
| 241 | + clj(1,j)*dfield(ij,1) & |
---|
| 242 | + clj(2,j)*dfield(ij,2) |
---|
| 243 | END DO |
---|
| 244 | CASE(3) |
---|
| 245 | !$OMP DO SCHEDULE(static) |
---|
| 246 | DO ij=1,sz |
---|
| 247 | field(ij) = field(ij) & |
---|
| 248 | + clj(1,j)*dfield(ij,1) & |
---|
| 249 | + clj(2,j)*dfield(ij,2) & |
---|
| 250 | + clj(3,j)*dfield(ij,3) |
---|
| 251 | |
---|
| 252 | END DO |
---|
| 253 | CASE(4) |
---|
| 254 | !$OMP DO SCHEDULE(static) |
---|
| 255 | DO ij=1,sz |
---|
| 256 | field(ij) = field(ij) & |
---|
| 257 | + clj(1,j)*dfield(ij,1) & |
---|
| 258 | + clj(2,j)*dfield(ij,2) & |
---|
| 259 | + clj(3,j)*dfield(ij,3) & |
---|
| 260 | + clj(4,j)*dfield(ij,4) |
---|
| 261 | END DO |
---|
| 262 | CASE(5) |
---|
| 263 | !$OMP DO SCHEDULE(static) |
---|
| 264 | DO ij=1,sz |
---|
| 265 | field(ij) = field(ij) & |
---|
| 266 | + clj(1,j)*dfield(ij,1) & |
---|
| 267 | + clj(2,j)*dfield(ij,2) & |
---|
| 268 | + clj(3,j)*dfield(ij,3) & |
---|
| 269 | + clj(4,j)*dfield(ij,4) & |
---|
| 270 | + clj(5,j)*dfield(ij,5) |
---|
| 271 | END DO |
---|
| 272 | END SELECT |
---|
[658] | 273 | CALL exit_trace() |
---|
[642] | 274 | END SUBROUTINE update |
---|
| 275 | |
---|
| 276 | !---------------------------------------------- XIOS ---------------------------------------- |
---|
[638] | 277 | |
---|
| 278 | #ifdef CPP_USING_XIOS |
---|
[642] | 279 | |
---|
[638] | 280 | SUBROUTINE setup_xios() BINDC(setup_xios) |
---|
| 281 | ! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine |
---|
| 282 | ! This is the case when calling from Python after importing mpi4py |
---|
| 283 | INTEGER :: ierr, mpi_threading_mode |
---|
| 284 | |
---|
| 285 | PRINT *, 'Initialize XIOS and obtain our communicator' |
---|
| 286 | CALL xios_initialize("icosagcm",return_comm=comm_icosa) |
---|
| 287 | |
---|
| 288 | PRINT *, 'Initialize our XIOS context' |
---|
| 289 | |
---|
| 290 | CALL xios_context_initialize("icosagcm",comm_icosa) |
---|
| 291 | CALL xios_get_handle("icosagcm",ctx_hdl) |
---|
| 292 | CALL xios_set_current_context(ctx_hdl) |
---|
| 293 | ! CALL xios_set_axis_attr("lev",n_glo=llm ,value=lev_value) ; |
---|
| 294 | ! CALL xios_set_axis_attr("levp1",n_glo=llm+1 ,value=lev_valuep1) ; |
---|
| 295 | ! CALL xios_set_axis_attr("nq",n_glo=nqtot, value=nq_value) ; |
---|
| 296 | ! CALL xios_set_domaingroup_attr("i",ni_glo=ncell_tot, ibegin=displ, ni=ncell) |
---|
| 297 | ! CALL xios_set_domaingroup_attr("i", data_dim=1, type='unstructured' , nvertex=6, i_index=ind_glo) |
---|
| 298 | ! CALL xios_set_domaingroup_attr("i",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat) |
---|
| 299 | ! CALL xios_set_domain_attr("u",ni_glo=ncell_tot, ibegin=displ, ni=ncell) |
---|
| 300 | ! CALL xios_set_domain_attr("u", data_dim=1, type='unstructured' , nvertex=2, i_index=ind_glo) |
---|
| 301 | ! CALL xios_set_domain_attr("u",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat) |
---|
| 302 | ! CALL xios_set_domain_attr("v",ni_glo=ncell_tot, ibegin=displ, ni=ncell) |
---|
| 303 | ! CALL xios_set_domain_attr("v", data_dim=1, type='unstructured' , nvertex=3) |
---|
| 304 | ! CALL xios_set_domain_attr("v",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat) |
---|
| 305 | ! CALL xios_set_timestep(dtime) |
---|
| 306 | ! CALL xios_set_fieldgroup_attr("standard_output", freq_op=itau_out*xios_timestep, freq_offset=(itau_out-1)*xios_timestep) |
---|
| 307 | ! CALL xios_close_context_definition() |
---|
| 308 | |
---|
| 309 | ! PRINT *, 'Read data' |
---|
| 310 | ! CALL xios_recv_field(name,field) |
---|
| 311 | |
---|
| 312 | ! PRINT *,'Write data' |
---|
| 313 | ! CALL xios_send_field(name,field_tmp) |
---|
| 314 | |
---|
| 315 | ! PRINT *, 'Flush to disk, clean up and die' |
---|
| 316 | ! CALL xios_context_finalize |
---|
| 317 | END SUBROUTINE setup_xios |
---|
| 318 | ! |
---|
| 319 | SUBROUTINE call_xios_set_timestep(dt) BINDC(xios_set_timestep) |
---|
| 320 | DBL, VALUE :: dt |
---|
| 321 | TYPE(xios_duration) :: dtime |
---|
| 322 | dtime%second=dt |
---|
| 323 | CALL xios_set_timestep(dtime) |
---|
| 324 | END SUBROUTINE call_xios_set_timestep |
---|
| 325 | |
---|
| 326 | SUBROUTINE call_xios_update_calendar(step) BINDC(xios_update_calendar) |
---|
| 327 | INDEX, value :: step |
---|
| 328 | CALL xios_update_calendar(step) |
---|
| 329 | END SUBROUTINE call_xios_update_calendar |
---|
| 330 | |
---|
| 331 | #endif |
---|
| 332 | |
---|
| 333 | END MODULE timestep_unstructured_mod |
---|