Changeset 2021 for LMDZ5/trunk/libf/dyn3d_common
- Timestamp:
- Apr 25, 2014, 12:20:14 PM (10 years ago)
- Location:
- LMDZ5/trunk/libf/dyn3d_common
- Files:
-
- 2 moved
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/dyn3d_common/exner_hyb_m.F90
r1992 r2021 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 -
LMDZ5/trunk/libf/dyn3d_common/exner_milieu_m.F90
r1992 r2021 1 ! 2 ! $Id $ 3 ! 4 SUBROUTINE exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ) 5 c 6 c Auteurs : F. Forget , Y. Wanherdrick 7 c P.Le Van , Fr. Hourdin . 8 c .......... 9 c 10 c .... ngrid, ps,p sont des argum.d'entree au sous-prog ... 11 c .... beta, pks,pk,pkf sont des argum.de sortie au sous-prog ... 12 c 13 c ************************************************************************ 14 c Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 15 c couches . Pk(l) sera calcule aux milieux des couches l ,entre les 16 c pressions p(l) et p(l+1) ,definis aux interfaces des llm couches . 17 c ************************************************************************ 18 c .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont 19 c la pression et la fonction d'Exner au sol . 20 c 21 c WARNING : CECI est une version speciale de exner_hyb originale 22 c Utilise dans la version martienne pour pouvoir 23 c tourner avec des coordonnees verticales complexe 24 c => Il ne verifie PAS la condition la proportionalite en 25 c energie totale/ interne / potentielle (F.Forget 2001) 26 c ( voir note de Fr.Hourdin ) , 27 c 28 IMPLICIT NONE 29 c 30 #include "dimensions.h" 31 #include "paramet.h" 32 #include "comconst.h" 33 #include "comgeom.h" 34 #include "comvert.h" 35 #include "serre.h" 1 module exner_milieu_m 36 2 37 INTEGER ngrid 38 REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm) 39 REAL ps(ngrid),pks(ngrid), beta(ngrid,llm) 3 IMPLICIT NONE 40 4 41 c .... variables locales ...5 contains 42 6 43 INTEGER l, ij 44 REAL dum1 7 SUBROUTINE exner_milieu ( ngrid, ps, p, pks, pk, pkf ) 8 ! 9 ! Auteurs : F. Forget , Y. Wanherdrick 10 ! P.Le Van , Fr. Hourdin . 11 ! .......... 12 ! 13 ! .... ngrid, ps,p sont des argum.d'entree au sous-prog ... 14 ! .... pks,pk,pkf sont des argum.de sortie au sous-prog ... 15 ! 16 ! ************************************************************************ 17 ! Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 18 ! couches . Pk(l) sera calcule aux milieux des couches l ,entre les 19 ! pressions p(l) et p(l+1) ,definis aux interfaces des llm couches . 20 ! ************************************************************************ 21 ! .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont 22 ! la pression et la fonction d'Exner au sol . 23 ! 24 ! WARNING : CECI est une version speciale de exner_hyb originale 25 ! Utilise dans la version martienne pour pouvoir 26 ! tourner avec des coordonnees verticales complexe 27 ! => Il ne verifie PAS la condition la proportionalite en 28 ! energie totale/ interne / potentielle (F.Forget 2001) 29 ! ( voir note de Fr.Hourdin ) , 30 ! 31 ! 32 include "dimensions.h" 33 include "paramet.h" 34 include "comconst.h" 35 include "comgeom.h" 36 include "comvert.h" 37 include "serre.h" 45 38 46 REAL ppn(iim),pps(iim) 47 REAL xpn, xps 48 REAL SSUM 49 EXTERNAL SSUM 50 logical,save :: firstcall=.true. 51 character(len=*),parameter :: modname="exner_milieu" 39 INTEGER ngrid 40 REAL p(ngrid,llmp1),pk(ngrid,llm) 41 real, optional:: pkf(ngrid,llm) 42 REAL ps(ngrid),pks(ngrid) 52 43 53 ! Sanity check 54 if (firstcall) then 55 ! sanity checks for Shallow Water case (1 vertical layer) 56 if (llm.eq.1) then 44 ! .... variables locales ... 45 46 INTEGER l, ij 47 REAL dum1 48 49 logical,save :: firstcall=.true. 50 character(len=*),parameter :: modname="exner_milieu" 51 52 ! Sanity check 53 if (firstcall) then 54 ! sanity checks for Shallow Water case (1 vertical layer) 55 if (llm.eq.1) then 57 56 if (kappa.ne.1) then 58 call abort_gcm(modname,59 &"kappa!=1 , but running in Shallow Water mode!!",42)57 call abort_gcm(modname, & 58 "kappa!=1 , but running in Shallow Water mode!!",42) 60 59 endif 61 60 if (cpp.ne.r) then 62 call abort_gcm(modname,63 &"cpp!=r , but running in Shallow Water mode!!",42)61 call abort_gcm(modname, & 62 "cpp!=r , but running in Shallow Water mode!!",42) 64 63 endif 65 64 endif ! of if (llm.eq.1) 66 65 67 68 66 firstcall=.false. 67 endif ! of if (firstcall) 69 68 70 !!!! Specific behaviour for Shallow Water (1 vertical layer) case:71 72 73 74 75 76 pks(ij) = (cpp/preff) * ps(ij) 69 ! Specific behaviour for Shallow Water (1 vertical layer) case: 70 if (llm.eq.1) then 71 72 ! Compute pks(:),pk(:),pkf(:) 73 74 DO ij = 1, ngrid 75 pks(ij) = (cpp/preff) * ps(ij) 77 76 pk(ij,1) = .5*pks(ij) 78 ENDDO 79 80 CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) 81 CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 82 83 ! our work is done, exit routine 84 return 77 ENDDO 85 78 86 endif ! of if (llm.eq.1) 79 if (present(pkf)) then 80 pkf = pk 81 CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 82 end if 87 83 88 !!!! General case: 84 ! our work is done, exit routine 85 return 86 endif ! of if (llm.eq.1) 89 87 90 c ------------- 91 c Calcul de pks 92 c ------------- 93 94 DO ij = 1, ngrid 95 pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 96 ENDDO 88 ! General case: 97 89 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 90 ! ------------- 91 ! Calcul de pks 92 ! ------------- 104 93 105 DO ij = 1, iip1 106 pks( ij ) = xpn 107 pks( ij+ip1jm ) = xps 108 ENDDO 109 c 110 c 111 c .... Calcul de pk pour la couche l 112 c -------------------------------------------- 113 c 114 dum1 = cpp * (2*preff)**(-kappa) 115 DO l = 1, llm-1 116 DO ij = 1, ngrid 117 pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa 118 ENDDO 119 ENDDO 94 DO ij = 1, ngrid 95 pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 96 ENDDO 120 97 121 c .... Calcul de pk pour la couche l = llm .. 122 c (on met la meme distance (en log pression) entre Pk(llm) 123 c et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2) 98 ! .... Calcul de pk pour la couche l 99 ! -------------------------------------------- 100 ! 101 dum1 = cpp * (2*preff)**(-kappa) 102 DO l = 1, llm-1 103 DO ij = 1, ngrid 104 pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa 105 ENDDO 106 ENDDO 124 107 125 DO ij = 1, ngrid126 pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)127 ENDDO108 ! .... Calcul de pk pour la couche l = llm .. 109 ! (on met la meme distance (en log pression) entre Pk(llm) 110 ! et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2) 128 111 112 DO ij = 1, ngrid 113 pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2) 114 ENDDO 129 115 130 c calcul de pkf 131 c ------------- 132 CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) 133 CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 134 135 c EST-CE UTILE ?? : calcul de beta 136 c -------------------------------- 137 DO l = 2, llm 138 DO ij = 1, ngrid 139 beta(ij,l) = pk(ij,l) / pk(ij,l-1) 140 ENDDO 141 ENDDO 116 if (present(pkf)) then 117 ! calcul de pkf 118 pkf = pk 119 CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 120 end if 142 121 143 RETURN 144 END 122 END SUBROUTINE exner_milieu 123 124 end module exner_milieu_m
Note: See TracChangeset
for help on using the changeset viewer.