| 1 | ! |
|---|
| 2 | ! $Id: inidissip.F90 1502 2011-03-21 16:07:54Z jghattas $ |
|---|
| 3 | ! |
|---|
| 4 | SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , & |
|---|
| 5 | tetagdiv,tetagrot,tetatemp ) |
|---|
| 6 | !======================================================================= |
|---|
| 7 | ! initialisation de la dissipation horizontale |
|---|
| 8 | !======================================================================= |
|---|
| 9 | !----------------------------------------------------------------------- |
|---|
| 10 | ! declarations: |
|---|
| 11 | ! ------------- |
|---|
| 12 | |
|---|
| 13 | USE control_mod, only : dissip_period,iperiod |
|---|
| 14 | |
|---|
| 15 | IMPLICIT NONE |
|---|
| 16 | include "dimensions.h" |
|---|
| 17 | include "paramet.h" |
|---|
| 18 | include "comdissipn.h" |
|---|
| 19 | include "comconst.h" |
|---|
| 20 | include "comvert.h" |
|---|
| 21 | include "logic.h" |
|---|
| 22 | include "iniprint.h" |
|---|
| 23 | |
|---|
| 24 | LOGICAL,INTENT(in) :: lstardis |
|---|
| 25 | INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh |
|---|
| 26 | REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp |
|---|
| 27 | |
|---|
| 28 | ! Local variables: |
|---|
| 29 | REAL fact,zvert(llm),zz |
|---|
| 30 | REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1) |
|---|
| 31 | real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm) |
|---|
| 32 | REAL ullm,vllm,umin,vmin,zhmin,zhmax |
|---|
| 33 | REAL zllm |
|---|
| 34 | |
|---|
| 35 | INTEGER l,ij,idum,ii |
|---|
| 36 | REAL tetamin |
|---|
| 37 | REAL Pup |
|---|
| 38 | character (len=80) :: abort_message |
|---|
| 39 | |
|---|
| 40 | REAL ran1 |
|---|
| 41 | |
|---|
| 42 | |
|---|
| 43 | !----------------------------------------------------------------------- |
|---|
| 44 | ! |
|---|
| 45 | ! calcul des valeurs propres des operateurs par methode iterrative: |
|---|
| 46 | ! ----------------------------------------------------------------- |
|---|
| 47 | |
|---|
| 48 | crot = -1. |
|---|
| 49 | cdivu = -1. |
|---|
| 50 | cdivh = -1. |
|---|
| 51 | |
|---|
| 52 | ! calcul de la valeur propre de divgrad: |
|---|
| 53 | ! -------------------------------------- |
|---|
| 54 | idum = 0 |
|---|
| 55 | DO l = 1, llm |
|---|
| 56 | DO ij = 1, ip1jmp1 |
|---|
| 57 | deltap(ij,l) = 1. |
|---|
| 58 | ENDDO |
|---|
| 59 | ENDDO |
|---|
| 60 | |
|---|
| 61 | idum = -1 |
|---|
| 62 | zh(1) = RAN1(idum)-.5 |
|---|
| 63 | idum = 0 |
|---|
| 64 | DO ij = 2, ip1jmp1 |
|---|
| 65 | zh(ij) = RAN1(idum) -.5 |
|---|
| 66 | ENDDO |
|---|
| 67 | |
|---|
| 68 | CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) |
|---|
| 69 | |
|---|
| 70 | CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) |
|---|
| 71 | |
|---|
| 72 | IF ( zhmin .GE. zhmax ) THEN |
|---|
| 73 | write(lunout,*)' Inidissip zh min max ',zhmin,zhmax |
|---|
| 74 | abort_message='probleme generateur alleatoire dans inidissip' |
|---|
| 75 | call abort_gcm('inidissip',abort_message,1) |
|---|
| 76 | ENDIF |
|---|
| 77 | |
|---|
| 78 | zllm = ABS( zhmax ) |
|---|
| 79 | DO l = 1,50 |
|---|
| 80 | IF(lstardis) THEN |
|---|
| 81 | CALL divgrad2(1,zh,deltap,niterh,divgra) |
|---|
| 82 | ELSE |
|---|
| 83 | CALL divgrad (1,zh,niterh,divgra) |
|---|
| 84 | ENDIF |
|---|
| 85 | |
|---|
| 86 | zllm = ABS(maxval(divgra)) |
|---|
| 87 | zh = divgra / zllm |
|---|
| 88 | ENDDO |
|---|
| 89 | |
|---|
| 90 | IF(lstardis) THEN |
|---|
| 91 | cdivh = 1./ zllm |
|---|
| 92 | ELSE |
|---|
| 93 | cdivh = zllm ** ( -1./niterh ) |
|---|
| 94 | ENDIF |
|---|
| 95 | |
|---|
| 96 | ! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) |
|---|
| 97 | ! ----------------------------------------------------------------- |
|---|
| 98 | write(lunout,*)'inidissip: calcul des valeurs propres' |
|---|
| 99 | |
|---|
| 100 | DO ii = 1, 2 |
|---|
| 101 | ! |
|---|
| 102 | DO ij = 1, ip1jmp1 |
|---|
| 103 | zu(ij) = RAN1(idum) -.5 |
|---|
| 104 | ENDDO |
|---|
| 105 | CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) |
|---|
| 106 | DO ij = 1, ip1jm |
|---|
| 107 | zv(ij) = RAN1(idum) -.5 |
|---|
| 108 | ENDDO |
|---|
| 109 | CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) |
|---|
| 110 | |
|---|
| 111 | CALL minmax(iip1*jjp1,zu,umin,ullm ) |
|---|
| 112 | CALL minmax(iip1*jjm, zv,vmin,vllm ) |
|---|
| 113 | |
|---|
| 114 | ullm = ABS ( ullm ) |
|---|
| 115 | vllm = ABS ( vllm ) |
|---|
| 116 | |
|---|
| 117 | DO l = 1, 50 |
|---|
| 118 | IF(ii.EQ.1) THEN |
|---|
| 119 | !cccc CALL covcont( 1,zu,zv,zu,zv ) |
|---|
| 120 | IF(lstardis) THEN |
|---|
| 121 | CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy ) |
|---|
| 122 | ELSE |
|---|
| 123 | CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy ) |
|---|
| 124 | ENDIF |
|---|
| 125 | ELSE |
|---|
| 126 | IF(lstardis) THEN |
|---|
| 127 | CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy ) |
|---|
| 128 | ELSE |
|---|
| 129 | CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy ) |
|---|
| 130 | ENDIF |
|---|
| 131 | ENDIF |
|---|
| 132 | |
|---|
| 133 | zllm = max(abs(maxval(gx)), abs(maxval(gy))) |
|---|
| 134 | zu = gx / zllm |
|---|
| 135 | zv = gy / zllm |
|---|
| 136 | end DO |
|---|
| 137 | |
|---|
| 138 | IF ( ii.EQ.1 ) THEN |
|---|
| 139 | IF(lstardis) THEN |
|---|
| 140 | cdivu = 1./zllm |
|---|
| 141 | ELSE |
|---|
| 142 | cdivu = zllm **( -1./nitergdiv ) |
|---|
| 143 | ENDIF |
|---|
| 144 | ELSE |
|---|
| 145 | IF(lstardis) THEN |
|---|
| 146 | crot = 1./ zllm |
|---|
| 147 | ELSE |
|---|
| 148 | crot = zllm **( -1./nitergrot ) |
|---|
| 149 | ENDIF |
|---|
| 150 | ENDIF |
|---|
| 151 | |
|---|
| 152 | end DO |
|---|
| 153 | |
|---|
| 154 | ! petit test pour les operateurs non star: |
|---|
| 155 | ! ---------------------------------------- |
|---|
| 156 | |
|---|
| 157 | ! IF(.NOT.lstardis) THEN |
|---|
| 158 | fact = rad*24./REAL(jjm) |
|---|
| 159 | fact = fact*fact |
|---|
| 160 | write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu |
|---|
| 161 | write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot |
|---|
| 162 | write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh |
|---|
| 163 | ! ENDIF |
|---|
| 164 | |
|---|
| 165 | !----------------------------------------------------------------------- |
|---|
| 166 | ! variation verticale du coefficient de dissipation: |
|---|
| 167 | ! -------------------------------------------------- |
|---|
| 168 | |
|---|
| 169 | ! First step: going from 1 to dissip_fac_mid (in gcm.def) |
|---|
| 170 | !============ |
|---|
| 171 | DO l=1,llm |
|---|
| 172 | zz = 1. - preff/presnivs(l) |
|---|
| 173 | zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) |
|---|
| 174 | ENDDO |
|---|
| 175 | |
|---|
| 176 | write(lunout,*) 'Dissipation : ' |
|---|
| 177 | write(lunout,*) 'Multiplication de la dissipation en altitude :' |
|---|
| 178 | write(lunout,*) ' dissip_fac_mid =', dissip_fac_mid |
|---|
| 179 | |
|---|
| 180 | ! Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def) |
|---|
| 181 | !========================== |
|---|
| 182 | ! Utilisation de la fonction tangente hyperbolique pour augmenter |
|---|
| 183 | ! arbitrairement la dissipation et donc la stabilite du modele a |
|---|
| 184 | ! partir d'une certaine altitude. |
|---|
| 185 | |
|---|
| 186 | ! Le facteur multiplicatif de basse atmosphere etant deja pris |
|---|
| 187 | ! en compte, il faut diviser le facteur multiplicatif de haute |
|---|
| 188 | ! atmosphere par celui-ci. |
|---|
| 189 | |
|---|
| 190 | if (ok_strato) then |
|---|
| 191 | |
|---|
| 192 | Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta) |
|---|
| 193 | do l=1,llm |
|---|
| 194 | zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) & |
|---|
| 195 | *(1-(0.5*(1+tanh(-6./dissip_deltaz* & |
|---|
| 196 | (-dissip_hdelta*log(presnivs(l)/Pup)) )))) )) |
|---|
| 197 | enddo |
|---|
| 198 | |
|---|
| 199 | write(*,*) ' dissip_fac_up =', dissip_fac_up |
|---|
| 200 | write(*,*) 'Transition mid /up: Pupstart,delta =', & |
|---|
| 201 | dissip_pupstart,'Pa', dissip_deltaz , '(km)' |
|---|
| 202 | |
|---|
| 203 | endif |
|---|
| 204 | |
|---|
| 205 | |
|---|
| 206 | write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale' |
|---|
| 207 | |
|---|
| 208 | tetamin = 1.e+6 |
|---|
| 209 | |
|---|
| 210 | DO l=1,llm |
|---|
| 211 | tetaudiv(l) = zvert(l)/tetagdiv |
|---|
| 212 | tetaurot(l) = zvert(l)/tetagrot |
|---|
| 213 | tetah(l) = zvert(l)/tetatemp |
|---|
| 214 | |
|---|
| 215 | IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) |
|---|
| 216 | IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) |
|---|
| 217 | IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) |
|---|
| 218 | ENDDO |
|---|
| 219 | |
|---|
| 220 | ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def |
|---|
| 221 | IF (dissip_period == 0) THEN |
|---|
| 222 | dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod |
|---|
| 223 | write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period |
|---|
| 224 | dissip_period = MAX(iperiod,dissip_period) |
|---|
| 225 | END IF |
|---|
| 226 | |
|---|
| 227 | dtdiss = dissip_period * dtvr |
|---|
| 228 | write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr |
|---|
| 229 | |
|---|
| 230 | DO l = 1,llm |
|---|
| 231 | write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), & |
|---|
| 232 | dtdiss*tetah(l) |
|---|
| 233 | ENDDO |
|---|
| 234 | |
|---|
| 235 | END SUBROUTINE inidissip |
|---|
