Changeset 1472 for LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F90
- Timestamp:
- Dec 23, 2010, 5:38:42 PM (14 years ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F90
r1463 r1472 1 !2 1 ! $Id$ 3 !4 SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)5 2 6 c Auteur : P. Le Van . 7 c 8 IMPLICIT NONE 3 SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) 9 4 10 #include "dimensions.h" 11 #include "paramet.h" 12 #include "iniprint.h" 13 #include "logic.h" 14 c 15 c======================================================================= 16 c 17 c 18 c s = sigma ** kappa : coordonnee verticale 19 c dsig(l) : epaisseur de la couche l ds la coord. s 20 c sig(l) : sigma a l'interface des couches l et l-1 21 c ds(l) : distance entre les couches l et l-1 en coord.s 22 c 23 c======================================================================= 24 c 25 REAL pa,preff 26 REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1) 27 REAL presnivs(llm) 28 c 29 c declarations: 30 c ------------- 31 c 32 REAL sig(llm+1),dsig(llm) 33 real zzz(1:llm+1) 34 real dzz(1:llm) 35 real zk,zkm1,dzk1,dzk2,k0,k1 36 c 37 INTEGER l 38 REAL snorm,dsigmin 39 REAL alpha,beta,gama,delta,deltaz,h 40 INTEGER np,ierr 41 REAL pi,x 5 ! Auteur : P. Le Van 42 6 43 REAL SSUM 44 c 45 c----------------------------------------------------------------------- 46 c 47 pi=2.*ASIN(1.) 7 IMPLICIT NONE 48 8 49 OPEN(99,file='sigma.def',status='old',form='formatted', 50 s iostat=ierr) 9 include "dimensions.h" 10 include "paramet.h" 11 include "iniprint.h" 12 include "logic.h" 51 13 52 c----------------------------------------------------------------------- 53 c cas 1 on lit les options dans sigma.def: 54 c ---------------------------------------- 14 ! s = sigma ** kappa : coordonnee verticale 15 ! dsig(l) : epaisseur de la couche l ds la coord. s 16 ! sig(l) : sigma a l'interface des couches l et l-1 17 ! ds(l) : distance entre les couches l et l-1 en coord.s 55 18 56 IF (ierr.eq.0) THEN 19 REAL pa, preff 20 REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1) 21 REAL presnivs(llm) 57 22 58 READ(99,*) h ! hauteur d'echelle 8. 59 READ(99,*) deltaz ! epaiseur de la premiere couche 0.04 60 READ(99,*) beta ! facteur d'acroissement en haut 1.3 61 READ(99,*) k0 ! nombre de couches dans la transition surf 62 READ(99,*) k1 ! nombre de couches dans la transition haute 63 CLOSE(99) 64 alpha=deltaz/(llm*h) 65 write(lunout,*)'h,alpha,k0,k1,beta' 23 REAL sig(llm+1), dsig(llm) 24 real zk, zkm1, dzk1, dzk2, k0, k1 66 25 67 c read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2 26 INTEGER l 27 REAL dsigmin 28 REAL alpha, beta, deltaz, h 29 INTEGER iostat 30 REAL pi, x 68 31 69 alpha=deltaz/tanh(1./k0)*2. 70 zkm1=0. 71 sig(1)=1. 72 do l=1,llm 73 sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) 74 + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) 32 !----------------------------------------------------------------------- 33 34 pi = 2 * ASIN(1.) 35 36 OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat) 37 38 IF (iostat == 0) THEN 39 ! cas 1 on lit les options dans sigma.def: 40 READ(99, *) h ! hauteur d'echelle 8. 41 READ(99, *) deltaz ! epaiseur de la premiere couche 0.04 42 READ(99, *) beta ! facteur d'acroissement en haut 1.3 43 READ(99, *) k0 ! nombre de couches dans la transition surf 44 READ(99, *) k1 ! nombre de couches dans la transition haute 45 CLOSE(99) 46 alpha=deltaz/(llm*h) 47 write(lunout, *)'h, alpha, k0, k1, beta' 48 49 alpha=deltaz/tanh(1./k0)*2. 50 zkm1=0. 51 sig(1)=1. 52 do l=1, llm 53 sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) & 54 *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) 75 55 zk=-h*log(sig(l+1)) 76 56 77 57 dzk1=alpha*tanh(l/k0) 78 58 dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) 79 write(lunout, *)l,sig(l+1),zk,zk-zkm1,dzk1,dzk259 write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 80 60 zkm1=zk 81 61 enddo 82 62 83 63 sig(llm+1)=0. 84 64 85 c 86 DO 2 l = 1, llm 87 dsig(l) = sig(l)-sig(l+1) 88 2 CONTINUE 89 c 65 DO l = 1, llm 66 dsig(l) = sig(l)-sig(l+1) 67 end DO 68 ELSE 69 if (ok_strato) then 70 if (llm==39) then 71 dsigmin=0.3 72 else if (llm==50) then 73 dsigmin=1. 74 else 75 WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster' 76 dsigmin=1. 77 endif 78 WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin 79 endif 90 80 91 ELSE 92 c----------------------------------------------------------------------- 93 c cas 2 ancienne discretisation (LMD5...): 94 c ---------------------------------------- 81 DO l = 1, llm 82 x = pi * (l - 0.5) / (llm + 1) 95 83 96 WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale' 84 IF (ok_strato) THEN 85 dsig(l) =(dsigmin + 7. * SIN(x)**2) & 86 *(1. - tanh(2 * x / pi - 1.))**2 / 4. 87 ELSE 88 dsig(l) = 1.0 + 7.0 * SIN(x)**2 89 ENDIF 90 ENDDO 91 dsig = dsig / sum(dsig) 92 sig(llm+1) = 0. 93 DO l = llm, 1, -1 94 sig(l) = sig(l+1) + dsig(l) 95 ENDDO 96 ENDIF 97 97 98 if (ok_strato) then 99 if (llm==39) then 100 dsigmin=0.3 101 else if (llm==50) then 102 dsigmin=1. 103 else 104 WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster' 105 dsigmin=1. 106 endif 107 WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin 108 endif 98 DO l=1, llm 99 nivsigs(l) = REAL(l) 100 ENDDO 109 101 110 h=7. 111 snorm = 0. 112 DO l = 1, llm 113 x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1) 102 DO l=1, llmp1 103 nivsig(l)= REAL(l) 104 ENDDO 114 105 115 IF (ok_strato) THEN 116 dsig(l) =(dsigmin + 7.0 * SIN(x)**2) 117 & *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 118 ELSE 119 dsig(l) = 1.0 + 7.0 * SIN(x)**2 120 ENDIF 106 ! .... Calculs de ap(l) et de bp(l) .... 107 ! ..... pa et preff sont lus sur les fichiers start par lectba ..... 121 108 122 snorm = snorm + dsig(l) 123 ENDDO 124 snorm = 1./snorm 125 DO l = 1, llm 126 dsig(l) = dsig(l)*snorm 127 ENDDO 128 sig(llm+1) = 0. 129 DO l = llm, 1, -1 130 sig(l) = sig(l+1) + dsig(l) 131 ENDDO 109 bp(llmp1) = 0. 132 110 133 ENDIF 111 DO l = 1, llm 112 bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) 113 ap(l) = pa * ( sig(l) - bp(l) ) 114 ENDDO 134 115 116 bp(1)=1. 117 ap(1)=0. 135 118 136 DO l=1,llm 137 nivsigs(l) = REAL(l) 138 ENDDO 119 ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) 139 120 140 DO l=1,llmp1 141 nivsig(l)= REAL(l) 142 ENDDO 121 write(lunout, *)' BP ' 122 write(lunout, *) bp 123 write(lunout, *)' AP ' 124 write(lunout, *) ap 143 125 144 c 145 c .... Calculs de ap(l) et de bp(l) .... 146 c ......................................... 147 c 148 c 149 c ..... pa et preff sont lus sur les fichiers start par lectba ..... 150 c 126 write(lunout, *) 'Niveaux de pressions approximatifs aux centres des' 127 write(lunout, *)'couches calcules pour une pression de surface =', preff 128 write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de' 129 write(lunout, *)'8km' 130 DO l = 1, llm 131 dpres(l) = bp(l) - bp(l+1) 132 presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) 133 write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & 134 log(preff/presnivs(l))*8. & 135 , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ & 136 max(ap(l+1)+bp(l+1)*preff, 1.e-10)) 137 ENDDO 151 138 152 bp(llmp1) = 0. 139 write(lunout, *) 'PRESNIVS ' 140 write(lunout, *) presnivs 153 141 154 DO l = 1, llm 155 cc 156 ccc ap(l) = 0. 157 ccc bp(l) = sig(l) 158 159 bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) 160 ap(l) = pa * ( sig(l) - bp(l) ) 161 c 162 ENDDO 163 164 bp(1)=1. 165 ap(1)=0. 166 167 ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) 168 169 write(lunout,*)' BP ' 170 write(lunout,*) bp 171 write(lunout,*)' AP ' 172 write(lunout,*) ap 173 174 write(lunout,*) 175 .'Niveaux de pressions approximatifs aux centres des' 176 write(lunout,*)'couches calcules pour une pression de surface =', 177 . preff 178 write(lunout,*) 179 . 'et altitudes equivalentes pour une hauteur d echelle de' 180 write(lunout,*)'8km' 181 DO l = 1, llm 182 dpres(l) = bp(l) - bp(l+1) 183 presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) 184 write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),' Z ~ ', 185 . log(preff/presnivs(l))*8. 186 . ,' DZ ~ ',8.*log((ap(l)+bp(l)*preff)/ 187 . max(ap(l+1)+bp(l+1)*preff,1.e-10)) 188 ENDDO 189 190 write(lunout,*)' PRESNIVS ' 191 write(lunout,*)presnivs 192 193 RETURN 194 END 142 END SUBROUTINE disvert
Note: See TracChangeset
for help on using the changeset viewer.