[12] | 1 | MODULE caldyn_gcm_mod |
---|
[19] | 2 | USE icosa |
---|
[151] | 3 | USE transfert_mod |
---|
[362] | 4 | USE caldyn_kernels_base_mod |
---|
[361] | 5 | USE caldyn_kernels_mod |
---|
[350] | 6 | IMPLICIT NONE |
---|
[132] | 7 | PRIVATE |
---|
| 8 | |
---|
[362] | 9 | PUBLIC init_caldyn, caldyn_BC, caldyn |
---|
| 10 | |
---|
[12] | 11 | CONTAINS |
---|
[15] | 12 | |
---|
[98] | 13 | SUBROUTINE init_caldyn |
---|
[50] | 14 | USE icosa |
---|
[354] | 15 | USE observable_mod |
---|
[131] | 16 | USE mpipara |
---|
[295] | 17 | USE omp_para |
---|
[50] | 18 | IMPLICIT NONE |
---|
[122] | 19 | CHARACTER(len=255) :: def |
---|
[179] | 20 | INTEGER :: ind |
---|
| 21 | REAL(rstd),POINTER :: planetvel(:) |
---|
[122] | 22 | |
---|
[195] | 23 | def='energy' |
---|
[125] | 24 | CALL getin('caldyn_conserv',def) |
---|
| 25 | SELECT CASE(TRIM(def)) |
---|
| 26 | CASE('energy') |
---|
[132] | 27 | caldyn_conserv=energy |
---|
[126] | 28 | CASE('enstrophy') |
---|
[132] | 29 | caldyn_conserv=enstrophy |
---|
[125] | 30 | CASE DEFAULT |
---|
[159] | 31 | IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', & |
---|
| 32 | TRIM(def),'> options are <energy>, <enstrophy>' |
---|
[125] | 33 | STOP |
---|
| 34 | END SELECT |
---|
[295] | 35 | IF (is_master) PRINT *, 'caldyn_conserv=',def |
---|
[125] | 36 | |
---|
[17] | 37 | CALL allocate_caldyn |
---|
[179] | 38 | |
---|
| 39 | DO ind=1,ndomain |
---|
[202] | 40 | IF (.NOT. assigned_domain(ind)) CYCLE |
---|
[179] | 41 | CALL swap_dimensions(ind) |
---|
| 42 | CALL swap_geometry(ind) |
---|
| 43 | planetvel=f_planetvel(ind) |
---|
| 44 | CALL compute_planetvel(planetvel) |
---|
| 45 | END DO |
---|
| 46 | |
---|
[15] | 47 | END SUBROUTINE init_caldyn |
---|
[179] | 48 | |
---|
[12] | 49 | SUBROUTINE allocate_caldyn |
---|
[19] | 50 | USE icosa |
---|
[12] | 51 | IMPLICIT NONE |
---|
| 52 | |
---|
[151] | 53 | CALL allocate_field(f_out_u,field_u,type_real,llm) |
---|
| 54 | CALL allocate_field(f_qu,field_u,type_real,llm) |
---|
[162] | 55 | CALL allocate_field(f_qv,field_z,type_real,llm) |
---|
[159] | 56 | CALL allocate_field(f_pk, field_t,type_real,llm, name='pk') |
---|
| 57 | CALL allocate_field(f_wwuu, field_u,type_real,llm+1,name='wwuu') |
---|
[179] | 58 | CALL allocate_field(f_planetvel, field_u,type_real, name='planetvel') ! planetary velocity at r=a |
---|
[151] | 59 | |
---|
[12] | 60 | END SUBROUTINE allocate_caldyn |
---|
[157] | 61 | |
---|
[350] | 62 | SUBROUTINE caldyn_BC(f_phis, f_geopot, f_wflux) |
---|
[157] | 63 | USE icosa |
---|
| 64 | USE mpipara |
---|
| 65 | USE omp_para |
---|
| 66 | TYPE(t_field),POINTER :: f_phis(:) |
---|
[350] | 67 | TYPE(t_field),POINTER :: f_geopot(:) |
---|
[157] | 68 | TYPE(t_field),POINTER :: f_wflux(:) |
---|
| 69 | REAL(rstd),POINTER :: phis(:) |
---|
| 70 | REAL(rstd),POINTER :: wflux(:,:) |
---|
| 71 | REAL(rstd),POINTER :: geopot(:,:) |
---|
| 72 | REAL(rstd),POINTER :: wwuu(:,:) |
---|
| 73 | |
---|
| 74 | INTEGER :: ind,i,j,ij,l |
---|
| 75 | |
---|
[295] | 76 | IF (is_omp_first_level) THEN |
---|
[157] | 77 | DO ind=1,ndomain |
---|
[186] | 78 | IF (.NOT. assigned_domain(ind)) CYCLE |
---|
[157] | 79 | CALL swap_dimensions(ind) |
---|
| 80 | CALL swap_geometry(ind) |
---|
| 81 | geopot=f_geopot(ind) |
---|
| 82 | phis=f_phis(ind) |
---|
| 83 | wflux=f_wflux(ind) |
---|
| 84 | wwuu=f_wwuu(ind) |
---|
| 85 | |
---|
[174] | 86 | DO ij=ij_begin_ext,ij_end_ext |
---|
| 87 | ! lower BCs : geopot=phis, wflux=0, wwuu=0 |
---|
| 88 | geopot(ij,1) = phis(ij) |
---|
| 89 | wflux(ij,1) = 0. |
---|
| 90 | wwuu(ij+u_right,1)=0 |
---|
| 91 | wwuu(ij+u_lup,1)=0 |
---|
| 92 | wwuu(ij+u_ldown,1)=0 |
---|
| 93 | ! top BCs : wflux=0, wwuu=0 |
---|
| 94 | wflux(ij,llm+1) = 0. |
---|
| 95 | wwuu(ij+u_right,llm+1)=0 |
---|
| 96 | wwuu(ij+u_lup,llm+1)=0 |
---|
| 97 | wwuu(ij+u_ldown,llm+1)=0 |
---|
[157] | 98 | ENDDO |
---|
| 99 | END DO |
---|
| 100 | ENDIF |
---|
| 101 | |
---|
[295] | 102 | !$OMP BARRIER |
---|
[157] | 103 | END SUBROUTINE caldyn_BC |
---|
[56] | 104 | |
---|
[159] | 105 | SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & |
---|
[350] | 106 | f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du) |
---|
[126] | 107 | USE icosa |
---|
[354] | 108 | USE observable_mod |
---|
[167] | 109 | USE disvert_mod, ONLY : caldyn_eta, eta_mass |
---|
[126] | 110 | USE vorticity_mod |
---|
| 111 | USE kinetic_mod |
---|
| 112 | USE theta2theta_rhodz_mod |
---|
[191] | 113 | USE wind_mod |
---|
[131] | 114 | USE mpipara |
---|
[145] | 115 | USE trace |
---|
[151] | 116 | USE omp_para |
---|
[171] | 117 | USE output_field_mod |
---|
[295] | 118 | USE checksum_mod |
---|
[126] | 119 | IMPLICIT NONE |
---|
[129] | 120 | LOGICAL,INTENT(IN) :: write_out |
---|
[126] | 121 | TYPE(t_field),POINTER :: f_phis(:) |
---|
[12] | 122 | TYPE(t_field),POINTER :: f_ps(:) |
---|
[159] | 123 | TYPE(t_field),POINTER :: f_mass(:) |
---|
[126] | 124 | TYPE(t_field),POINTER :: f_theta_rhodz(:) |
---|
| 125 | TYPE(t_field),POINTER :: f_u(:) |
---|
| 126 | TYPE(t_field),POINTER :: f_q(:) |
---|
[350] | 127 | TYPE(t_field),POINTER :: f_geopot(:) |
---|
[134] | 128 | TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:) |
---|
[360] | 129 | TYPE(t_field) :: f_dps(:) |
---|
| 130 | TYPE(t_field) :: f_dmass(:) |
---|
| 131 | TYPE(t_field) :: f_dtheta_rhodz(:) |
---|
| 132 | TYPE(t_field) :: f_du(:) |
---|
[12] | 133 | |
---|
[159] | 134 | REAL(rstd),POINTER :: ps(:), dps(:) |
---|
[387] | 135 | REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:) |
---|
[134] | 136 | REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:) |
---|
[159] | 137 | REAL(rstd),POINTER :: qu(:,:) |
---|
[162] | 138 | REAL(rstd),POINTER :: qv(:,:) |
---|
[151] | 139 | |
---|
| 140 | ! temporary shared variable |
---|
| 141 | REAL(rstd),POINTER :: theta(:,:) |
---|
[157] | 142 | REAL(rstd),POINTER :: pk(:,:) |
---|
[156] | 143 | REAL(rstd),POINTER :: geopot(:,:) |
---|
[162] | 144 | REAL(rstd),POINTER :: convm(:,:) |
---|
[151] | 145 | REAL(rstd),POINTER :: wwuu(:,:) |
---|
[157] | 146 | |
---|
| 147 | INTEGER :: ind |
---|
[151] | 148 | LOGICAL,SAVE :: first=.TRUE. |
---|
| 149 | !$OMP THREADPRIVATE(first) |
---|
[12] | 150 | |
---|
[162] | 151 | IF (first) THEN |
---|
[151] | 152 | first=.FALSE. |
---|
[162] | 153 | IF(caldyn_eta==eta_mass) THEN |
---|
| 154 | CALL init_message(f_ps,req_i1,req_ps) |
---|
[190] | 155 | ELSE |
---|
[162] | 156 | CALL init_message(f_mass,req_i1,req_mass) |
---|
| 157 | END IF |
---|
[151] | 158 | CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz) |
---|
| 159 | CALL init_message(f_u,req_e1_vect,req_u) |
---|
| 160 | CALL init_message(f_qu,req_e1_scal,req_qu) |
---|
[361] | 161 | ! Overlapping com/compute (deactivated) : MPI messages need to be sent at first call to caldyn |
---|
| 162 | ! This is needed only once : the next ones will be sent by timeloop |
---|
[186] | 163 | ! IF(caldyn_eta==eta_mass) THEN |
---|
| 164 | ! CALL send_message(f_ps,req_ps) |
---|
| 165 | ! CALL wait_message(req_ps) |
---|
| 166 | ! ELSE |
---|
| 167 | ! CALL send_message(f_mass,req_mass) |
---|
| 168 | ! CALL wait_message(req_mass) |
---|
| 169 | ! END IF |
---|
| 170 | ENDIF |
---|
| 171 | |
---|
| 172 | CALL trace_start("caldyn") |
---|
| 173 | |
---|
[361] | 174 | IF(caldyn_eta==eta_mass) THEN |
---|
| 175 | CALL send_message(f_ps,req_ps) ! COM00 |
---|
| 176 | ELSE |
---|
| 177 | CALL send_message(f_mass,req_mass) ! COM00 |
---|
| 178 | END IF |
---|
| 179 | |
---|
| 180 | CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01 |
---|
| 181 | CALL send_message(f_u,req_u) ! COM02 |
---|
[126] | 182 | |
---|
| 183 | SELECT CASE(caldyn_conserv) |
---|
[132] | 184 | CASE(energy) ! energy-conserving |
---|
[128] | 185 | DO ind=1,ndomain |
---|
[186] | 186 | IF (.NOT. assigned_domain(ind)) CYCLE |
---|
[128] | 187 | CALL swap_dimensions(ind) |
---|
| 188 | CALL swap_geometry(ind) |
---|
| 189 | ps=f_ps(ind) |
---|
[157] | 190 | u=f_u(ind) |
---|
| 191 | theta_rhodz = f_theta_rhodz(ind) |
---|
[159] | 192 | mass=f_mass(ind) |
---|
[157] | 193 | theta = f_theta(ind) |
---|
[128] | 194 | qu=f_qu(ind) |
---|
[162] | 195 | qv=f_qv(ind) |
---|
[387] | 196 | CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv) ! COM00 COM01 COM02 |
---|
[128] | 197 | ENDDO |
---|
[347] | 198 | ! CALL checksum(f_mass) |
---|
| 199 | ! CALL checksum(f_theta) |
---|
[128] | 200 | |
---|
[361] | 201 | CALL send_message(f_qu,req_qu) ! COM03 wait_message in caldyn_horiz |
---|
[327] | 202 | ! CALL wait_message(req_qu) |
---|
[128] | 203 | |
---|
| 204 | DO ind=1,ndomain |
---|
[186] | 205 | IF (.NOT. assigned_domain(ind)) CYCLE |
---|
[128] | 206 | CALL swap_dimensions(ind) |
---|
| 207 | CALL swap_geometry(ind) |
---|
| 208 | ps=f_ps(ind) |
---|
[157] | 209 | u=f_u(ind) |
---|
[159] | 210 | mass=f_mass(ind) |
---|
[157] | 211 | theta = f_theta(ind) |
---|
[128] | 212 | qu=f_qu(ind) |
---|
[151] | 213 | pk = f_pk(ind) |
---|
[156] | 214 | geopot = f_geopot(ind) |
---|
[183] | 215 | CALL compute_geopot(ps,mass,theta, pk,geopot) |
---|
[157] | 216 | hflux=f_hflux(ind) |
---|
[162] | 217 | convm = f_dmass(ind) |
---|
[157] | 218 | dtheta_rhodz=f_dtheta_rhodz(ind) |
---|
| 219 | du=f_du(ind) |
---|
[387] | 220 | CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du) |
---|
[162] | 221 | IF(caldyn_eta==eta_mass) THEN |
---|
| 222 | wflux=f_wflux(ind) |
---|
| 223 | wwuu=f_wwuu(ind) |
---|
| 224 | dps=f_dps(ind) |
---|
[387] | 225 | CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz(:,:,1), du) |
---|
[162] | 226 | END IF |
---|
[128] | 227 | ENDDO |
---|
[347] | 228 | |
---|
| 229 | ! CALL checksum(f_geopot) |
---|
| 230 | ! CALL checksum(f_dmass) |
---|
| 231 | ! CALL checksum(f_pk) |
---|
| 232 | ! CALL checksum(f_pk) |
---|
[128] | 233 | |
---|
[132] | 234 | CASE(enstrophy) ! enstrophy-conserving |
---|
[126] | 235 | DO ind=1,ndomain |
---|
[186] | 236 | IF (.NOT. assigned_domain(ind)) CYCLE |
---|
[126] | 237 | CALL swap_dimensions(ind) |
---|
| 238 | CALL swap_geometry(ind) |
---|
[157] | 239 | ps=f_ps(ind) |
---|
| 240 | u=f_u(ind) |
---|
| 241 | theta_rhodz=f_theta_rhodz(ind) |
---|
[159] | 242 | mass=f_mass(ind) |
---|
[157] | 243 | theta = f_theta(ind) |
---|
| 244 | qu=f_qu(ind) |
---|
[182] | 245 | qv=f_qv(ind) |
---|
[387] | 246 | CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv) |
---|
[157] | 247 | pk = f_pk(ind) |
---|
[156] | 248 | geopot = f_geopot(ind) |
---|
[183] | 249 | CALL compute_geopot(ps,mass,theta, pk,geopot) |
---|
[134] | 250 | hflux=f_hflux(ind) |
---|
[162] | 251 | convm = f_dmass(ind) |
---|
[126] | 252 | dtheta_rhodz=f_dtheta_rhodz(ind) |
---|
| 253 | du=f_du(ind) |
---|
[387] | 254 | CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du) |
---|
[162] | 255 | IF(caldyn_eta==eta_mass) THEN |
---|
| 256 | wflux=f_wflux(ind) |
---|
| 257 | wwuu=f_wwuu(ind) |
---|
| 258 | dps=f_dps(ind) |
---|
| 259 | CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du) |
---|
| 260 | END IF |
---|
[126] | 261 | ENDDO |
---|
| 262 | |
---|
| 263 | CASE DEFAULT |
---|
| 264 | STOP |
---|
| 265 | END SELECT |
---|
[12] | 266 | |
---|
[295] | 267 | !$OMP BARRIER |
---|
[126] | 268 | ! CALL check_mass_conservation(f_ps,f_dps) |
---|
[145] | 269 | CALL trace_end("caldyn") |
---|
[186] | 270 | !!$OMP BARRIER |
---|
[126] | 271 | |
---|
| 272 | END SUBROUTINE caldyn |
---|
[179] | 273 | |
---|
[126] | 274 | !-------------------------------- Diagnostics ---------------------------- |
---|
| 275 | |
---|
| 276 | SUBROUTINE check_mass_conservation(f_ps,f_dps) |
---|
| 277 | USE icosa |
---|
[131] | 278 | USE mpipara |
---|
[126] | 279 | IMPLICIT NONE |
---|
| 280 | TYPE(t_field),POINTER :: f_ps(:) |
---|
| 281 | TYPE(t_field),POINTER :: f_dps(:) |
---|
| 282 | REAL(rstd),POINTER :: ps(:) |
---|
| 283 | REAL(rstd),POINTER :: dps(:) |
---|
| 284 | REAL(rstd) :: mass_tot,dmass_tot |
---|
| 285 | INTEGER :: ind,i,j,ij |
---|
| 286 | |
---|
| 287 | mass_tot=0 |
---|
| 288 | dmass_tot=0 |
---|
| 289 | |
---|
| 290 | CALL transfert_request(f_dps,req_i1) |
---|
| 291 | CALL transfert_request(f_ps,req_i1) |
---|
| 292 | |
---|
| 293 | DO ind=1,ndomain |
---|
| 294 | CALL swap_dimensions(ind) |
---|
| 295 | CALL swap_geometry(ind) |
---|
| 296 | |
---|
| 297 | ps=f_ps(ind) |
---|
| 298 | dps=f_dps(ind) |
---|
| 299 | |
---|
| 300 | DO j=jj_begin,jj_end |
---|
| 301 | DO i=ii_begin,ii_end |
---|
| 302 | ij=(j-1)*iim+i |
---|
| 303 | IF (domain(ind)%own(i,j)) THEN |
---|
| 304 | mass_tot=mass_tot+ps(ij)*Ai(ij)/g |
---|
| 305 | dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g |
---|
| 306 | ENDIF |
---|
| 307 | ENDDO |
---|
| 308 | ENDDO |
---|
| 309 | |
---|
| 310 | ENDDO |
---|
[131] | 311 | IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot," dmass_tot ",dmass_tot |
---|
[126] | 312 | |
---|
| 313 | END SUBROUTINE check_mass_conservation |
---|
[12] | 314 | |
---|
| 315 | END MODULE caldyn_gcm_mod |
---|