Changeset 2021 for LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc_m.F90
- Timestamp:
- Apr 25, 2014, 12:20:14 PM (11 years ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc_m.F90
r1992 r2021 1 c 2 c $Id$ 3 c 4 SUBROUTINE exner_hyb_loc(ngrid, ps, p,alpha,beta, pks,pk,pkf) 5 c 6 c Auteurs : P.Le Van , Fr. Hourdin . 7 c .......... 8 c 9 c .... ngrid, ps,p sont des argum.d'entree au sous-prog ... 10 c .... alpha,beta, pks,pk,pkf sont des argum.de sortie au sous-prog ... 11 c 12 c ************************************************************************ 13 c Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des 14 c couches . Pk(l) sera calcule aux milieux des couches l ,entre les 15 c pressions p(l) et p(l+1) ,definis aux interfaces des llm couches . 16 c ************************************************************************ 17 c .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont 18 c la pression et la fonction d'Exner au sol . 19 c 20 c -------- z 21 c A partir des relations ( 1 ) p*dz(pk) = kappa *pk*dz(p) et 22 c ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1) 23 c ( voir note de Fr.Hourdin ) , 24 c 25 c on determine successivement , du haut vers le bas des couches, les 26 c coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 27 c puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 28 c pk(ij,l) donne par la relation (2), pour l = 2 a l = llm . 29 c 30 c 31 USE parallel_lmdz 32 USE mod_filtreg_p 33 USE write_field_loc 34 IMPLICIT NONE 35 c 36 #include "dimensions.h" 37 #include "paramet.h" 38 #include "comconst.h" 39 #include "comgeom.h" 40 #include "comvert.h" 41 #include "serre.h" 1 module exner_hyb_loc_m 42 2 43 INTEGER ngrid 44 REAL p(ijb_u:ije_u,llmp1),pk(ijb_u:ije_u,llm) 45 REAL pkf(ijb_u:ije_u,llm) 46 REAL ps(ijb_u:ije_u),pks(ijb_u:ije_u) 47 REAL alpha(ijb_u:ije_u,llm),beta(ijb_u:ije_u,llm) 3 IMPLICIT NONE 48 4 49 c .... variables locales ...5 contains 50 6 51 INTEGER l, ij 52 REAL unpl2k,dellta 7 SUBROUTINE exner_hyb_loc(ngrid, ps, p, pks,pk,pkf) 53 8 54 REAL ppn(iim),pps(iim) 55 REAL xpn, xps 56 REAL SSUM 57 EXTERNAL SSUM 58 INTEGER ije,ijb,jje,jjb 59 logical,save :: firstcall=.true. 60 !$OMP THREADPRIVATE(firstcall) 61 character(len=*),parameter :: modname="exner_hyb_loc" 62 c 63 c$OMP BARRIER 9 ! Auteurs : P.Le Van , Fr. Hourdin . 10 ! .......... 11 ! 12 ! .... ngrid, ps,p sont des argum.d'entree au sous-prog ... 13 ! .... pks,pk,pkf sont des argum.de sortie au sous-prog ... 14 ! 15 ! ************************************************************************ 16 ! Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 17 ! couches . Pk(l) sera calcule aux milieux des couches l ,entre les 18 ! pressions p(l) et p(l+1) ,definis aux interfaces des llm couches . 19 ! ************************************************************************ 20 ! .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont 21 ! la pression et la fonction d'Exner au sol . 22 ! 23 ! -------- z 24 ! A partir des relations ( 1 ) p*dz(pk) = kappa *pk*dz(p) et 25 ! ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1) 26 ! ( voir note de Fr.Hourdin ) , 27 ! 28 ! on determine successivement , du haut vers le bas des couches, les 29 ! coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 30 ! puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 31 ! pk(ij,l) donne par la relation (2), pour l = 2 a l = llm . 32 ! 33 ! 34 USE parallel_lmdz 35 USE mod_filtreg_p 36 USE write_field_loc 37 ! 38 include "dimensions.h" 39 include "paramet.h" 40 include "comconst.h" 41 include "comgeom.h" 42 include "comvert.h" 43 include "serre.h" 64 44 65 ! Sanity check 66 if (firstcall) then 67 ! sanity checks for Shallow Water case (1 vertical layer) 68 if (llm.eq.1) then 45 INTEGER ngrid 46 REAL p(ijb_u:ije_u,llmp1),pk(ijb_u:ije_u,llm) 47 REAL, optional:: pkf(ijb_u:ije_u,llm) 48 REAL ps(ijb_u:ije_u),pks(ijb_u:ije_u) 49 REAL alpha(ijb_u:ije_u,llm),beta(ijb_u:ije_u,llm) 50 51 ! .... variables locales ... 52 53 INTEGER l, ij 54 REAL unpl2k,dellta 55 56 INTEGER ije,ijb,jje,jjb 57 logical,save :: firstcall=.true. 58 !$OMP THREADPRIVATE(firstcall) 59 character(len=*),parameter :: modname="exner_hyb_loc" 60 ! 61 !$OMP BARRIER 62 63 ! Sanity check 64 if (firstcall) then 65 ! sanity checks for Shallow Water case (1 vertical layer) 66 if (llm.eq.1) then 69 67 if (kappa.ne.1) then 70 call abort_gcm(modname,71 &"kappa!=1 , but running in Shallow Water mode!!",42)68 call abort_gcm(modname, & 69 "kappa!=1 , but running in Shallow Water mode!!",42) 72 70 endif 73 71 if (cpp.ne.r) then 74 call abort_gcm(modname,75 &"cpp!=r , but running in Shallow Water mode!!",42)72 call abort_gcm(modname, & 73 "cpp!=r , but running in Shallow Water mode!!",42) 76 74 endif 77 75 endif ! of if (llm.eq.1) 78 76 79 80 77 firstcall=.false. 78 endif ! of if (firstcall) 81 79 82 c$OMP BARRIER80 !$OMP BARRIER 83 81 84 ! Specific behaviour for Shallow Water (1 vertical layer) case 85 86 87 88 89 90 !$OMP DO SCHEDULE(STATIC)91 92 pks(ij) =(cpp/preff)*ps(ij)82 ! Specific behaviour for Shallow Water (1 vertical layer) case: 83 if (llm.eq.1) then 84 85 ! Compute pks(:),pk(:),pkf(:) 86 ijb=ij_begin 87 ije=ij_end 88 !$OMP DO SCHEDULE(STATIC) 89 DO ij=ijb, ije 90 pks(ij) = (cpp/preff) * ps(ij) 93 91 pk(ij,1) = .5*pks(ij) 94 pkf(ij,1)=pk(ij,1)95 96 !$OMP ENDDO92 if (present(pkf)) pkf(ij,1)=pk(ij,1) 93 ENDDO 94 !$OMP ENDDO 97 95 98 !$OMP MASTER 99 if (pole_nord) then 100 DO ij = 1, iim 101 ppn(ij) = aire( ij ) * pks( ij ) 102 ENDDO 103 xpn = SSUM(iim,ppn,1) /apoln 104 105 DO ij = 1, iip1 106 pks( ij ) = xpn 107 pk(ij,1) = .5*pks(ij) 108 pkf(ij,1)=pk(ij,1) 109 ENDDO 110 endif 111 112 if (pole_sud) then 113 DO ij = 1, iim 114 pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm ) 115 ENDDO 116 xps = SSUM(iim,pps,1) /apols 117 118 DO ij = 1, iip1 119 pks( ij+ip1jm ) = xps 120 pk(ij+ip1jm,1)=.5*pks(ij+ip1jm) 121 pkf(ij+ip1jm,1)=pk(ij+ip1jm,1) 122 ENDDO 123 endif 124 !$OMP END MASTER 125 !$OMP BARRIER 126 jjb=jj_begin 127 jje=jj_end 128 CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, 129 & 2, 1, .TRUE., 1 ) 96 !$OMP BARRIER 97 if (present(pkf)) then 98 jjb=jj_begin 99 jje=jj_end 100 CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, & 101 2, 1, .TRUE., 1 ) 102 end if 130 103 131 132 133 104 ! our work is done, exit routine 105 return 106 endif ! of if (llm.eq.1) 134 107 108 ! General case: 135 109 136 unpl2k = 1.+ 2.* kappa 137 c 138 ijb=ij_begin 139 ije=ij_end 110 unpl2k = 1.+ 2.* kappa 140 111 141 c$OMP DO SCHEDULE(STATIC) 142 DO ij = ijb, ije 143 pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 144 ENDDO 145 c$OMP ENDDO 146 c Synchro OPENMP ici 112 ! ------------- 113 ! Calcul de pks 114 ! ------------- 147 115 148 c$OMP MASTER 149 if (pole_nord) then 150 DO ij = 1, iim 151 ppn(ij) = aire( ij ) * pks( ij ) 152 ENDDO 153 xpn = SSUM(iim,ppn,1) /apoln 154 155 DO ij = 1, iip1 156 pks( ij ) = xpn 157 ENDDO 158 endif 159 160 if (pole_sud) then 161 DO ij = 1, iim 162 pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm ) 163 ENDDO 164 xps = SSUM(iim,pps,1) /apols 165 166 DO ij = 1, iip1 167 pks( ij+ip1jm ) = xps 168 ENDDO 169 endif 170 c$OMP END MASTER 171 c$OMP BARRIER 172 c 173 c 174 c .... Calcul des coeff. alpha et beta pour la couche l = llm .. 175 c 176 c$OMP DO SCHEDULE(STATIC) 177 DO ij = ijb,ije 116 ijb=ij_begin 117 ije=ij_end 118 119 !$OMP DO SCHEDULE(STATIC) 120 DO ij = ijb, ije 121 pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 122 ENDDO 123 !$OMP ENDDO 124 ! Synchro OPENMP ici 125 126 !$OMP BARRIER 127 ! 128 ! 129 ! .... Calcul des coeff. alpha et beta pour la couche l = llm .. 130 ! 131 !$OMP DO SCHEDULE(STATIC) 132 DO ij = ijb,ije 178 133 alpha(ij,llm) = 0. 179 134 beta (ij,llm) = 1./ unpl2k 180 ENDDO 181 c$OMP ENDDO NOWAIT 182 c 183 c ... Calcul des coeff. alpha et beta pour l = llm-1 a l = 2 ... 184 c 185 DO l = llm -1 , 2 , -1 186 c 187 c$OMP DO SCHEDULE(STATIC) 188 DO ij = ijb, ije 189 dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k ) 190 alpha(ij,l) = - p(ij,l+1) / dellta * alpha(ij,l+1) 191 beta (ij,l) = p(ij,l ) / dellta 192 ENDDO 193 c$OMP ENDDO NOWAIT 194 c 195 ENDDO 135 ENDDO 136 !$OMP ENDDO NOWAIT 137 ! 138 ! ... Calcul des coeff. alpha et beta pour l = llm-1 a l = 2 ... 139 ! 140 DO l = llm -1 , 2 , -1 141 ! 142 !$OMP DO SCHEDULE(STATIC) 143 DO ij = ijb, ije 144 dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k ) 145 alpha(ij,l) = - p(ij,l+1) / dellta * alpha(ij,l+1) 146 beta (ij,l) = p(ij,l ) / dellta 147 ENDDO 148 !$OMP ENDDO NOWAIT 149 ENDDO 196 150 197 c 198 c *********************************************************************** 199 c ..... Calcul de pk pour la couche 1 , pres du sol .... 200 c 201 c$OMP DO SCHEDULE(STATIC) 202 DO ij = ijb, ije 203 pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) ) / 204 * ( p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) ) 205 ENDDO 206 c$OMP ENDDO NOWAIT 207 c 208 c ..... Calcul de pk(ij,l) , pour l = 2 a l = llm ........ 209 c 210 DO l = 2, llm 211 c$OMP DO SCHEDULE(STATIC) 212 DO ij = ijb, ije 213 pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1) 214 ENDDO 215 c$OMP ENDDO NOWAIT 216 ENDDO 217 c 218 c 219 c CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) 220 DO l = 1, llm 221 c$OMP DO SCHEDULE(STATIC) 222 DO ij = ijb, ije 223 pkf(ij,l)=pk(ij,l) 224 ENDDO 225 c$OMP ENDDO NOWAIT 226 ENDDO 151 ! *********************************************************************** 152 ! ..... Calcul de pk pour la couche 1 , pres du sol .... 153 ! 154 !$OMP DO SCHEDULE(STATIC) 155 DO ij = ijb, ije 156 pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) ) / & 157 ( p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) ) 158 ENDDO 159 !$OMP ENDDO NOWAIT 160 ! 161 ! ..... Calcul de pk(ij,l) , pour l = 2 a l = llm ........ 162 ! 163 DO l = 2, llm 164 !$OMP DO SCHEDULE(STATIC) 165 DO ij = ijb, ije 166 pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1) 167 ENDDO 168 !$OMP ENDDO NOWAIT 169 ENDDO 227 170 228 c$OMP BARRIER 229 230 jjb=jj_begin 231 jje=jj_end 171 if (present(pkf)) then 172 ! calcul de pkf 173 174 DO l = 1, llm 175 !$OMP DO SCHEDULE(STATIC) 176 DO ij = ijb, ije 177 pkf(ij,l)=pk(ij,l) 178 ENDDO 179 !$OMP ENDDO NOWAIT 180 ENDDO 181 182 !$OMP BARRIER 183 184 jjb=jj_begin 185 jje=jj_end 232 186 #ifdef DEBUG_IO 233 call WriteField_u('pkf',pkf)187 call WriteField_u('pkf',pkf) 234 188 #endif 235 CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm,236 &2, 1, .TRUE., 1 )189 CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, & 190 2, 1, .TRUE., 1 ) 237 191 #ifdef DEBUG_IO 238 call WriteField_u('pkf',pkf) 239 #endif 192 call WriteField_u('pkf',pkf) 193 #endif 194 end if 240 195 241 RETURN 242 END 196 END SUBROUTINE exner_hyb_loc 197 198 end module exner_hyb_loc_m
Note: See TracChangeset
for help on using the changeset viewer.