Changeset 1302 for trunk/LMDZ.COMMON/libf/dyn3d_common/exner_hyb_m.F90
- Timestamp:
- Jun 26, 2014, 6:07:05 PM (10 years ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/dyn3d_common/exner_hyb_m.F90
r1300 r1302 1 ! 2 ! $Id $ 3 ! 4 SUBROUTINE exner_hyb ( 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 IMPLICIT NONE 32 c 33 #include "dimensions.h" 34 #include "paramet.h" 35 #include "comconst.h" 36 #include "comgeom.h" 37 #include "comvert.h" 38 #include "serre.h" 1 module exner_hyb_m 39 2 40 INTEGER ngrid 41 REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm) 42 REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm) 3 IMPLICIT NONE 43 4 44 c .... variables locales ...5 contains 45 6 46 INTEGER l, ij 47 REAL unpl2k,dellta 7 SUBROUTINE exner_hyb ( ngrid, ps, p, pks, pk, pkf ) 48 8 49 REAL ppn(iim),pps(iim) 50 REAL xpn, xps 51 REAL SSUM 52 c 53 logical,save :: firstcall=.true. 54 character(len=*),parameter :: modname="exner_hyb" 55 56 ! Sanity check 57 if (firstcall) then 58 ! sanity checks for Shallow Water case (1 vertical layer) 59 if (llm.eq.1) then 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 ! 35 include "dimensions.h" 36 include "paramet.h" 37 include "comconst.h" 38 include "comgeom.h" 39 include "comvert.h" 40 include "serre.h" 41 42 INTEGER ngrid 43 REAL p(ngrid,llmp1),pk(ngrid,llm) 44 real, optional:: pkf(ngrid,llm) 45 REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm) 46 47 ! .... variables locales ... 48 49 INTEGER l, ij 50 REAL unpl2k,dellta 51 52 logical,save :: firstcall=.true. 53 character(len=*),parameter :: modname="exner_hyb" 54 55 ! Sanity check 56 if (firstcall) then 57 ! sanity checks for Shallow Water case (1 vertical layer) 58 if (llm.eq.1) then 60 59 if (kappa.ne.1) then 61 call abort_gcm(modname,62 &"kappa!=1 , but running in Shallow Water mode!!",42)60 call abort_gcm(modname, & 61 "kappa!=1 , but running in Shallow Water mode!!",42) 63 62 endif 64 63 if (cpp.ne.r) then 65 call abort_gcm(modname,66 &"cpp!=r , but running in Shallow Water mode!!",42)64 call abort_gcm(modname, & 65 "cpp!=r , but running in Shallow Water mode!!",42) 67 66 endif 68 67 endif ! of if (llm.eq.1) 69 68 70 71 69 firstcall=.false. 70 endif ! of if (firstcall) 72 71 73 if (llm.eq.1) then 74 75 ! Compute pks(:),pk(:),pkf(:) 76 77 DO ij = 1, ngrid 78 pks(ij) = (cpp/preff) * ps(ij) 72 ! Specific behaviour for Shallow Water (1 vertical layer) case: 73 if (llm.eq.1) then 74 75 ! Compute pks(:),pk(:),pkf(:) 76 77 DO ij = 1, ngrid 78 pks(ij) = (cpp/preff) * ps(ij) 79 79 pk(ij,1) = .5*pks(ij) 80 ENDDO 81 82 CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) 83 CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 84 85 ! our work is done, exit routine 86 return 80 ENDDO 87 81 88 endif ! of if (llm.eq.1) 82 if (present(pkf)) then 83 pkf = pk 84 CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 85 end if 89 86 90 !!!! General case: 91 92 unpl2k = 1.+ 2.* kappa 93 c 94 DO ij = 1, ngrid 95 pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 96 ENDDO 87 ! our work is done, exit routine 88 return 89 endif ! of if (llm.eq.1) 97 90 98 DO ij = 1, iim 99 ppn(ij) = aire( ij ) * pks( ij ) 100 pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm ) 101 ENDDO 102 xpn = SSUM(iim,ppn,1) /apoln 103 xps = SSUM(iim,pps,1) /apols 91 ! General case: 104 92 105 DO ij = 1, iip1 106 pks( ij ) = xpn 107 pks( ij+ip1jm ) = xps 108 ENDDO 109 c 110 c 111 c .... Calcul des coeff. alpha et beta pour la couche l = llm .. 112 c 113 DO ij = 1, ngrid 93 unpl2k = 1.+ 2.* kappa 94 95 ! ------------- 96 ! Calcul de pks 97 ! ------------- 98 99 DO ij = 1, ngrid 100 pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 101 ENDDO 102 103 ! .... Calcul des coeff. alpha et beta pour la couche l = llm .. 104 ! 105 DO ij = 1, ngrid 114 106 alpha(ij,llm) = 0. 115 107 beta (ij,llm) = 1./ unpl2k 116 ENDDO 117 c 118 c ... Calcul des coeff. alpha et beta pour l = llm-1 a l = 2 ... 119 c 120 DO l = llm -1 , 2 , -1 121 c 122 DO ij = 1, ngrid 123 dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k ) 124 alpha(ij,l) = - p(ij,l+1) / dellta * alpha(ij,l+1) 125 beta (ij,l) = p(ij,l ) / dellta 126 ENDDO 127 c 128 ENDDO 129 c 130 c *********************************************************************** 131 c ..... Calcul de pk pour la couche 1 , pres du sol .... 132 c 133 DO ij = 1, ngrid 134 pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) ) / 135 * ( p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) ) 136 ENDDO 137 c 138 c ..... Calcul de pk(ij,l) , pour l = 2 a l = llm ........ 139 c 140 DO l = 2, llm 141 DO ij = 1, ngrid 142 pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1) 143 ENDDO 144 ENDDO 145 c 146 c 147 CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) 148 CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 149 108 ENDDO 109 ! 110 ! ... Calcul des coeff. alpha et beta pour l = llm-1 a l = 2 ... 111 ! 112 DO l = llm -1 , 2 , -1 113 ! 114 DO ij = 1, ngrid 115 dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k ) 116 alpha(ij,l) = - p(ij,l+1) / dellta * alpha(ij,l+1) 117 beta (ij,l) = p(ij,l ) / dellta 118 ENDDO 119 ENDDO 150 120 151 RETURN 152 END 121 ! *********************************************************************** 122 ! ..... Calcul de pk pour la couche 1 , pres du sol .... 123 ! 124 DO ij = 1, ngrid 125 pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) ) / & 126 ( p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) ) 127 ENDDO 128 ! 129 ! ..... Calcul de pk(ij,l) , pour l = 2 a l = llm ........ 130 ! 131 DO l = 2, llm 132 DO ij = 1, ngrid 133 pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1) 134 ENDDO 135 ENDDO 136 137 if (present(pkf)) then 138 ! calcul de pkf 139 pkf = pk 140 CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 141 end if 142 143 END SUBROUTINE exner_hyb 144 145 end module exner_hyb_m 146
Note: See TracChangeset
for help on using the changeset viewer.
