[1279] | 1 | ! $Id: disvert.F90 1645 2012-07-30 16:01:50Z idelkadi $ |
---|
[630] | 2 | |
---|
[1520] | 3 | SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight) |
---|
[630] | 4 | |
---|
[1472] | 5 | ! Auteur : P. Le Van |
---|
[630] | 6 | |
---|
[1620] | 7 | use new_unit_m, only: new_unit |
---|
| 8 | use ioipsl, only: getin |
---|
[1645] | 9 | use assert_m, only: assert |
---|
[1620] | 10 | |
---|
[1472] | 11 | IMPLICIT NONE |
---|
[630] | 12 | |
---|
[1472] | 13 | include "dimensions.h" |
---|
| 14 | include "paramet.h" |
---|
| 15 | include "iniprint.h" |
---|
| 16 | include "logic.h" |
---|
[630] | 17 | |
---|
[1472] | 18 | ! s = sigma ** kappa : coordonnee verticale |
---|
| 19 | ! dsig(l) : epaisseur de la couche l ds la coord. s |
---|
| 20 | ! sig(l) : sigma a l'interface des couches l et l-1 |
---|
| 21 | ! ds(l) : distance entre les couches l et l-1 en coord.s |
---|
[630] | 22 | |
---|
[1520] | 23 | real,intent(in) :: pa, preff |
---|
[1645] | 24 | real,intent(out) :: ap(llmp1) ! in Pa |
---|
| 25 | real, intent(out):: bp(llmp1) |
---|
[1520] | 26 | real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1) |
---|
| 27 | real,intent(out) :: presnivs(llm) |
---|
| 28 | real,intent(out) :: scaleheight |
---|
[630] | 29 | |
---|
[1472] | 30 | REAL sig(llm+1), dsig(llm) |
---|
| 31 | real zk, zkm1, dzk1, dzk2, k0, k1 |
---|
[630] | 32 | |
---|
[1620] | 33 | INTEGER l, unit |
---|
[1472] | 34 | REAL dsigmin |
---|
[1524] | 35 | REAL alpha, beta, deltaz |
---|
[1520] | 36 | REAL x |
---|
| 37 | character(len=*),parameter :: modname="disvert" |
---|
[630] | 38 | |
---|
[1620] | 39 | character(len=6):: vert_sampling |
---|
| 40 | ! (allowed values are "param", "tropo", "strato" and "read") |
---|
| 41 | |
---|
[1472] | 42 | !----------------------------------------------------------------------- |
---|
| 43 | |
---|
[1620] | 44 | print *, "Call sequence information: disvert" |
---|
| 45 | |
---|
[1520] | 46 | ! default scaleheight is 8km for earth |
---|
| 47 | scaleheight=8. |
---|
[1472] | 48 | |
---|
[1623] | 49 | vert_sampling = merge("strato", "tropo ", ok_strato) ! default value |
---|
[1620] | 50 | call getin('vert_sampling', vert_sampling) |
---|
| 51 | print *, 'vert_sampling = ' // vert_sampling |
---|
[1472] | 52 | |
---|
[1620] | 53 | select case (vert_sampling) |
---|
| 54 | case ("param") |
---|
| 55 | ! On lit les options dans sigma.def: |
---|
| 56 | OPEN(99, file='sigma.def', status='old', form='formatted') |
---|
[1520] | 57 | READ(99, *) scaleheight ! hauteur d'echelle 8. |
---|
[1472] | 58 | READ(99, *) deltaz ! epaiseur de la premiere couche 0.04 |
---|
| 59 | READ(99, *) beta ! facteur d'acroissement en haut 1.3 |
---|
| 60 | READ(99, *) k0 ! nombre de couches dans la transition surf |
---|
| 61 | READ(99, *) k1 ! nombre de couches dans la transition haute |
---|
| 62 | CLOSE(99) |
---|
[1520] | 63 | alpha=deltaz/(llm*scaleheight) |
---|
| 64 | write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', & |
---|
| 65 | scaleheight, alpha, k0, k1, beta |
---|
[1472] | 66 | |
---|
| 67 | alpha=deltaz/tanh(1./k0)*2. |
---|
| 68 | zkm1=0. |
---|
| 69 | sig(1)=1. |
---|
| 70 | do l=1, llm |
---|
[1520] | 71 | sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) & |
---|
| 72 | *exp(-alpha/scaleheight*tanh((llm-k1)/k0) & |
---|
| 73 | *beta**(l-(llm-k1))/log(beta)) |
---|
| 74 | zk=-scaleheight*log(sig(l+1)) |
---|
[630] | 75 | |
---|
| 76 | dzk1=alpha*tanh(l/k0) |
---|
| 77 | dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) |
---|
[1472] | 78 | write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 |
---|
[630] | 79 | zkm1=zk |
---|
[1472] | 80 | enddo |
---|
[630] | 81 | |
---|
[1472] | 82 | sig(llm+1)=0. |
---|
[1645] | 83 | |
---|
| 84 | bp(: llm) = EXP(1. - 1. / sig(: llm)**2) |
---|
| 85 | bp(llmp1) = 0. |
---|
| 86 | |
---|
| 87 | ap = pa * (sig - bp) |
---|
[1620] | 88 | case("tropo") |
---|
[1472] | 89 | DO l = 1, llm |
---|
[1620] | 90 | x = 2*asin(1.) * (l - 0.5) / (llm + 1) |
---|
| 91 | dsig(l) = 1.0 + 7.0 * SIN(x)**2 |
---|
| 92 | ENDDO |
---|
| 93 | dsig = dsig / sum(dsig) |
---|
| 94 | sig(llm+1) = 0. |
---|
| 95 | DO l = llm, 1, -1 |
---|
| 96 | sig(l) = sig(l+1) + dsig(l) |
---|
| 97 | ENDDO |
---|
[1645] | 98 | |
---|
| 99 | bp(1)=1. |
---|
| 100 | bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2) |
---|
| 101 | bp(llmp1) = 0. |
---|
| 102 | |
---|
| 103 | ap(1)=0. |
---|
| 104 | ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1)) |
---|
[1620] | 105 | case("strato") |
---|
| 106 | if (llm==39) then |
---|
| 107 | dsigmin=0.3 |
---|
| 108 | else if (llm==50) then |
---|
| 109 | dsigmin=1. |
---|
| 110 | else |
---|
| 111 | write(lunout,*) trim(modname), ' ATTENTION discretisation z a ajuster' |
---|
| 112 | dsigmin=1. |
---|
[1472] | 113 | endif |
---|
[1620] | 114 | WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin |
---|
[630] | 115 | |
---|
[1472] | 116 | DO l = 1, llm |
---|
[1480] | 117 | x = 2*asin(1.) * (l - 0.5) / (llm + 1) |
---|
[1620] | 118 | dsig(l) =(dsigmin + 7. * SIN(x)**2) & |
---|
| 119 | *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 |
---|
[1472] | 120 | ENDDO |
---|
| 121 | dsig = dsig / sum(dsig) |
---|
| 122 | sig(llm+1) = 0. |
---|
| 123 | DO l = llm, 1, -1 |
---|
| 124 | sig(l) = sig(l+1) + dsig(l) |
---|
| 125 | ENDDO |
---|
[1645] | 126 | |
---|
| 127 | bp(1)=1. |
---|
| 128 | bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2) |
---|
| 129 | bp(llmp1) = 0. |
---|
| 130 | |
---|
| 131 | ap(1)=0. |
---|
| 132 | ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1)) |
---|
[1620] | 133 | case("read") |
---|
[1645] | 134 | ! Read "ap" and "bp". First line is skipped (title line). "ap" |
---|
| 135 | ! should be in Pa. First couple of values should correspond to |
---|
| 136 | ! the surface, that is : "bp" should be in descending order. |
---|
[1620] | 137 | call new_unit(unit) |
---|
| 138 | open(unit, file="hybrid.txt", status="old", action="read", & |
---|
| 139 | position="rewind") |
---|
| 140 | read(unit, fmt=*) ! skip title line |
---|
| 141 | do l = 1, llm + 1 |
---|
[1645] | 142 | read(unit, fmt=*) ap(l), bp(l) |
---|
[1620] | 143 | end do |
---|
| 144 | close(unit) |
---|
[1645] | 145 | call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., & |
---|
| 146 | bp(llm + 1) == 0., "disvert: bad ap or bp values") |
---|
[1620] | 147 | case default |
---|
| 148 | call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1) |
---|
| 149 | END select |
---|
[1000] | 150 | |
---|
[1472] | 151 | DO l=1, llm |
---|
| 152 | nivsigs(l) = REAL(l) |
---|
| 153 | ENDDO |
---|
[630] | 154 | |
---|
[1472] | 155 | DO l=1, llmp1 |
---|
| 156 | nivsig(l)= REAL(l) |
---|
| 157 | ENDDO |
---|
[630] | 158 | |
---|
[1520] | 159 | write(lunout, *) trim(modname),': BP ' |
---|
[1472] | 160 | write(lunout, *) bp |
---|
[1520] | 161 | write(lunout, *) trim(modname),': AP ' |
---|
[1472] | 162 | write(lunout, *) ap |
---|
[630] | 163 | |
---|
[1472] | 164 | write(lunout, *) 'Niveaux de pressions approximatifs aux centres des' |
---|
| 165 | write(lunout, *)'couches calcules pour une pression de surface =', preff |
---|
[1520] | 166 | write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de ' |
---|
| 167 | write(lunout, *) scaleheight,' km' |
---|
[1472] | 168 | DO l = 1, llm |
---|
| 169 | dpres(l) = bp(l) - bp(l+1) |
---|
| 170 | presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) |
---|
| 171 | write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & |
---|
[1520] | 172 | log(preff/presnivs(l))*scaleheight & |
---|
| 173 | , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ & |
---|
[1472] | 174 | max(ap(l+1)+bp(l+1)*preff, 1.e-10)) |
---|
| 175 | ENDDO |
---|
[1279] | 176 | |
---|
[1520] | 177 | write(lunout, *) trim(modname),': PRESNIVS ' |
---|
[1472] | 178 | write(lunout, *) presnivs |
---|
[1279] | 179 | |
---|
[1472] | 180 | END SUBROUTINE disvert |
---|