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