[17] | 1 | MODULE caldyn_adv_mod |
---|
[19] | 2 | USE icosa |
---|
[17] | 3 | |
---|
[139] | 4 | IMPLICIT NONE |
---|
| 5 | PRIVATE |
---|
| 6 | PUBLIC :: init_caldyn, caldyn |
---|
| 7 | |
---|
[17] | 8 | CONTAINS |
---|
[139] | 9 | |
---|
[98] | 10 | SUBROUTINE init_caldyn |
---|
[17] | 11 | END SUBROUTINE init_caldyn |
---|
| 12 | |
---|
| 13 | SUBROUTINE check_mass_conservation(f_ps,f_dps) |
---|
[139] | 14 | USE icosa |
---|
[17] | 15 | TYPE(t_field),POINTER :: f_ps(:) |
---|
| 16 | TYPE(t_field),POINTER :: f_dps(:) |
---|
| 17 | REAL(rstd),POINTER :: ps(:) |
---|
| 18 | REAL(rstd),POINTER :: dps(:) |
---|
| 19 | REAL(rstd) :: mass_tot,dmass_tot |
---|
| 20 | INTEGER :: ind,i,j,ij |
---|
[139] | 21 | |
---|
[17] | 22 | mass_tot=0 |
---|
| 23 | dmass_tot=0 |
---|
[139] | 24 | |
---|
[17] | 25 | CALL transfert_request(f_dps,req_i1) |
---|
| 26 | CALL transfert_request(f_ps,req_i1) |
---|
| 27 | |
---|
| 28 | DO ind=1,ndomain |
---|
[139] | 29 | CALL swap_dimensions(ind) |
---|
| 30 | CALL swap_geometry(ind) |
---|
[17] | 31 | |
---|
[139] | 32 | ps=f_ps(ind) |
---|
| 33 | dps=f_dps(ind) |
---|
[17] | 34 | |
---|
[139] | 35 | DO j=jj_begin,jj_end |
---|
| 36 | DO i=ii_begin,ii_end |
---|
| 37 | ij=(j-1)*iim+i |
---|
| 38 | IF (domain(ind)%own(i,j)) THEN |
---|
| 39 | mass_tot=mass_tot+ps(ij)*Ai(ij)/g |
---|
| 40 | dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g |
---|
| 41 | ENDIF |
---|
| 42 | ENDDO |
---|
| 43 | ENDDO |
---|
| 44 | |
---|
[17] | 45 | ENDDO |
---|
| 46 | PRINT*, "mass_tot ", mass_tot," dmass_tot ",dmass_tot |
---|
| 47 | |
---|
| 48 | END SUBROUTINE check_mass_conservation |
---|
| 49 | |
---|
| 50 | |
---|
[139] | 51 | SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, & |
---|
| 52 | f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) |
---|
| 53 | USE icosa |
---|
| 54 | USE vorticity_mod |
---|
| 55 | USE kinetic_mod |
---|
| 56 | USE theta2theta_rhodz_mod |
---|
| 57 | IMPLICIT NONE |
---|
| 58 | LOGICAL,INTENT(IN) :: write_out |
---|
| 59 | TYPE(t_field),POINTER :: f_phis(:) |
---|
| 60 | TYPE(t_field),POINTER :: f_ps(:) |
---|
| 61 | TYPE(t_field),POINTER :: f_theta_rhodz(:) |
---|
| 62 | TYPE(t_field),POINTER :: f_u(:) |
---|
| 63 | TYPE(t_field),POINTER :: f_q(:) |
---|
| 64 | TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:) |
---|
| 65 | TYPE(t_field),POINTER :: f_dps(:) |
---|
| 66 | TYPE(t_field),POINTER :: f_dtheta_rhodz(:) |
---|
| 67 | TYPE(t_field),POINTER :: f_du(:) |
---|
| 68 | |
---|
| 69 | REAL(rstd),POINTER :: ps(:) |
---|
| 70 | REAL(rstd),POINTER :: u(:,:) |
---|
| 71 | REAL(rstd),POINTER :: dps(:) |
---|
| 72 | REAL(rstd),POINTER :: hflux(:,:), wflux(:,:) |
---|
| 73 | REAL(rstd),POINTER :: dtheta_rhodz(:,:), du(:,:) ! set to 0 |
---|
| 74 | INTEGER :: ind |
---|
| 75 | |
---|
[17] | 76 | CALL transfert_request(f_ps,req_i1) |
---|
[146] | 77 | CALL transfert_request(f_u,req_e1_vect) |
---|
[139] | 78 | |
---|
[17] | 79 | DO ind=1,ndomain |
---|
[139] | 80 | CALL swap_dimensions(ind) |
---|
| 81 | CALL swap_geometry(ind) |
---|
| 82 | ps=f_ps(ind) |
---|
| 83 | u=f_u(ind) |
---|
| 84 | dps=f_dps(ind) |
---|
| 85 | hflux=f_hflux(ind) |
---|
| 86 | wflux=f_wflux(ind) |
---|
| 87 | dtheta_rhodz=f_dtheta_rhodz(ind) |
---|
| 88 | du=f_du(ind) |
---|
| 89 | |
---|
| 90 | !$OMP PARALLEL DEFAULT(SHARED) |
---|
| 91 | CALL compute_caldyn(ps,u,hflux, wflux, dps, dtheta_rhodz, du) |
---|
| 92 | !$OMP END PARALLEL |
---|
[17] | 93 | ENDDO |
---|
| 94 | |
---|
[129] | 95 | IF (write_out) THEN |
---|
[139] | 96 | CALL writefield("ps",f_ps) |
---|
| 97 | CALL writefield("wflux",f_wflux) |
---|
[17] | 98 | ENDIF |
---|
[139] | 99 | ! CALL check_mass_conservation(f_ps,f_dps) |
---|
[17] | 100 | |
---|
| 101 | END SUBROUTINE caldyn |
---|
[139] | 102 | |
---|
| 103 | |
---|
| 104 | SUBROUTINE compute_caldyn(ps,u, hflux,wflux,dps, dtheta_rhodz,du) |
---|
| 105 | USE icosa |
---|
| 106 | USE disvert_mod |
---|
| 107 | IMPLICIT NONE |
---|
| 108 | REAL(rstd),INTENT(IN) :: ps(iim*jjm) |
---|
| 109 | REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm) |
---|
| 110 | REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) ! hflux in kg/s |
---|
| 111 | REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) |
---|
| 112 | REAL(rstd),INTENT(OUT) :: dps(iim*jjm) |
---|
| 113 | REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s) |
---|
| 114 | |
---|
| 115 | REAL(rstd),ALLOCATABLE :: rhodz(:,:) |
---|
| 116 | REAL(rstd),ALLOCATABLE :: divm(:,:) ! mass flux divergence |
---|
| 117 | |
---|
| 118 | INTEGER :: i,j,ij,l |
---|
[17] | 119 | LOGICAL,SAVE :: first=.TRUE. |
---|
[139] | 120 | |
---|
| 121 | ALLOCATE(rhodz(iim*jjm,llm)) |
---|
| 122 | ALLOCATE(divm(iim*jjm,llm)) ! mass flux divergence |
---|
| 123 | |
---|
| 124 | dtheta_rhodz(:,:)=0. |
---|
| 125 | du(:,:)=0. |
---|
| 126 | |
---|
[17] | 127 | !!! Compute mass |
---|
[139] | 128 | DO l = 1, llm |
---|
| 129 | DO j=jj_begin-1,jj_end+1 |
---|
| 130 | DO i=ii_begin-1,ii_end+1 |
---|
| 131 | ij=(j-1)*iim+i |
---|
| 132 | rhodz(ij,l) = (ap(l)-ap(l+1) + ps(ij)*(bp(l)-bp(l+1)) )/g |
---|
| 133 | ENDDO |
---|
[17] | 134 | ENDDO |
---|
[139] | 135 | ENDDO |
---|
[17] | 136 | |
---|
[139] | 137 | DO l = 1, llm |
---|
| 138 | !!! Mass fluxes |
---|
| 139 | DO j=jj_begin-1,jj_end+1 |
---|
| 140 | DO i=ii_begin-1,ii_end+1 |
---|
| 141 | ij=(j-1)*iim+i |
---|
| 142 | hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) |
---|
| 143 | hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) |
---|
| 144 | hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) |
---|
| 145 | ENDDO |
---|
| 146 | ENDDO |
---|
| 147 | !!! Horizontal divergence of fluxes |
---|
| 148 | DO j=jj_begin,jj_end |
---|
| 149 | DO i=ii_begin,ii_end |
---|
| 150 | ij=(j-1)*iim+i |
---|
| 151 | ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 |
---|
| 152 | divm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l) + & |
---|
| 153 | ne(ij,rup)*hflux(ij+u_rup,l) + & |
---|
| 154 | ne(ij,lup)*hflux(ij+u_lup,l) + & |
---|
| 155 | ne(ij,left)*hflux(ij+u_left,l) + & |
---|
| 156 | ne(ij,ldown)*hflux(ij+u_ldown,l) + & |
---|
| 157 | ne(ij,rdown)*hflux(ij+u_rdown,l)) |
---|
| 158 | ENDDO |
---|
| 159 | ENDDO |
---|
| 160 | ENDDO |
---|
[17] | 161 | |
---|
[139] | 162 | !!! cumulate mass flux divergence from top to bottom |
---|
| 163 | DO l = llm-1, 1, -1 |
---|
| 164 | !$OMP DO |
---|
| 165 | DO j=jj_begin,jj_end |
---|
| 166 | DO i=ii_begin,ii_end |
---|
| 167 | ij=(j-1)*iim+i |
---|
| 168 | divm(ij,l) = divm(ij,l) + divm(ij,l+1) |
---|
| 169 | ENDDO |
---|
| 170 | ENDDO |
---|
[17] | 171 | ENDDO |
---|
| 172 | |
---|
[139] | 173 | !!! Compute vertical mass flux |
---|
| 174 | DO l = 1,llm-1 |
---|
| 175 | DO j=jj_begin,jj_end |
---|
| 176 | DO i=ii_begin,ii_end |
---|
| 177 | ij=(j-1)*iim+i |
---|
| 178 | ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt |
---|
| 179 | ! => w>0 for upward transport |
---|
| 180 | wflux( ij, l+1 ) = divm( ij, l+1 ) - bp(l+1) * divm( ij, 1 ) |
---|
| 181 | ENDDO |
---|
| 182 | ENDDO |
---|
| 183 | ENDDO |
---|
[17] | 184 | |
---|
[139] | 185 | ! compute dps, set vertical mass flux at the surface to 0 |
---|
[17] | 186 | DO j=jj_begin,jj_end |
---|
[139] | 187 | DO i=ii_begin,ii_end |
---|
| 188 | ij=(j-1)*iim+i |
---|
| 189 | wflux(ij,1) = 0. |
---|
| 190 | ! dps/dt = -int(div flux)dz |
---|
| 191 | dps(ij)=-divm(ij,1) * g |
---|
[17] | 192 | ENDDO |
---|
| 193 | ENDDO |
---|
| 194 | |
---|
[139] | 195 | DEALLOCATE(rhodz) |
---|
| 196 | DEALLOCATE(divm) |
---|
| 197 | |
---|
[17] | 198 | END SUBROUTINE compute_caldyn |
---|
[139] | 199 | |
---|
[17] | 200 | END MODULE caldyn_adv_mod |
---|