[630] | 1 | ! |
---|
[1279] | 2 | ! $Id: inidissip.F90 1698 2012-12-19 17:14:07Z musat $ |
---|
[630] | 3 | ! |
---|
[1489] | 4 | SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , & |
---|
[1697] | 5 | tetagdiv,tetagrot,tetatemp, vert_prof_dissip) |
---|
[1489] | 6 | !======================================================================= |
---|
| 7 | ! initialisation de la dissipation horizontale |
---|
| 8 | !======================================================================= |
---|
| 9 | !----------------------------------------------------------------------- |
---|
| 10 | ! declarations: |
---|
| 11 | ! ------------- |
---|
[630] | 12 | |
---|
[1502] | 13 | USE control_mod, only : dissip_period,iperiod |
---|
[1403] | 14 | |
---|
[1489] | 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" |
---|
[1490] | 22 | include "iniprint.h" |
---|
[630] | 23 | |
---|
[1490] | 24 | LOGICAL,INTENT(in) :: lstardis |
---|
| 25 | INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh |
---|
| 26 | REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp |
---|
| 27 | |
---|
[1697] | 28 | integer, INTENT(in):: vert_prof_dissip |
---|
| 29 | ! Vertical profile of horizontal dissipation |
---|
| 30 | ! Allowed values: |
---|
[1698] | 31 | ! 0: rational fraction, function of pressure |
---|
[1697] | 32 | ! 1: tanh of altitude |
---|
| 33 | |
---|
[1490] | 34 | ! Local variables: |
---|
[1489] | 35 | REAL fact,zvert(llm),zz |
---|
[1611] | 36 | REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1) |
---|
| 37 | real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm) |
---|
[1489] | 38 | REAL ullm,vllm,umin,vmin,zhmin,zhmax |
---|
[1611] | 39 | REAL zllm |
---|
[630] | 40 | |
---|
[1489] | 41 | INTEGER l,ij,idum,ii |
---|
| 42 | REAL tetamin |
---|
| 43 | REAL pseudoz |
---|
[1490] | 44 | character (len=80) :: abort_message |
---|
[630] | 45 | |
---|
[1489] | 46 | REAL ran1 |
---|
[630] | 47 | |
---|
| 48 | |
---|
[1489] | 49 | !----------------------------------------------------------------------- |
---|
| 50 | ! |
---|
| 51 | ! calcul des valeurs propres des operateurs par methode iterrative: |
---|
| 52 | ! ----------------------------------------------------------------- |
---|
[630] | 53 | |
---|
[1489] | 54 | crot = -1. |
---|
| 55 | cdivu = -1. |
---|
| 56 | cdivh = -1. |
---|
[630] | 57 | |
---|
[1489] | 58 | ! calcul de la valeur propre de divgrad: |
---|
| 59 | ! -------------------------------------- |
---|
| 60 | idum = 0 |
---|
| 61 | DO l = 1, llm |
---|
| 62 | DO ij = 1, ip1jmp1 |
---|
[630] | 63 | deltap(ij,l) = 1. |
---|
[1489] | 64 | ENDDO |
---|
| 65 | ENDDO |
---|
[630] | 66 | |
---|
[1489] | 67 | idum = -1 |
---|
| 68 | zh(1) = RAN1(idum)-.5 |
---|
| 69 | idum = 0 |
---|
| 70 | DO ij = 2, ip1jmp1 |
---|
| 71 | zh(ij) = RAN1(idum) -.5 |
---|
| 72 | ENDDO |
---|
[630] | 73 | |
---|
[1489] | 74 | CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) |
---|
[630] | 75 | |
---|
[1489] | 76 | CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) |
---|
[630] | 77 | |
---|
[1489] | 78 | IF ( zhmin .GE. zhmax ) THEN |
---|
[1490] | 79 | write(lunout,*)' Inidissip zh min max ',zhmin,zhmax |
---|
| 80 | abort_message='probleme generateur alleatoire dans inidissip' |
---|
| 81 | call abort_gcm('inidissip',abort_message,1) |
---|
[1489] | 82 | ENDIF |
---|
[630] | 83 | |
---|
[1489] | 84 | zllm = ABS( zhmax ) |
---|
| 85 | DO l = 1,50 |
---|
| 86 | IF(lstardis) THEN |
---|
[1611] | 87 | CALL divgrad2(1,zh,deltap,niterh,divgra) |
---|
[1489] | 88 | ELSE |
---|
[1611] | 89 | CALL divgrad (1,zh,niterh,divgra) |
---|
[1489] | 90 | ENDIF |
---|
[630] | 91 | |
---|
[1611] | 92 | zllm = ABS(maxval(divgra)) |
---|
| 93 | zh = divgra / zllm |
---|
[1489] | 94 | ENDDO |
---|
[630] | 95 | |
---|
[1489] | 96 | IF(lstardis) THEN |
---|
| 97 | cdivh = 1./ zllm |
---|
| 98 | ELSE |
---|
| 99 | cdivh = zllm ** ( -1./niterh ) |
---|
| 100 | ENDIF |
---|
[630] | 101 | |
---|
[1489] | 102 | ! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) |
---|
| 103 | ! ----------------------------------------------------------------- |
---|
[1490] | 104 | write(lunout,*)'inidissip: calcul des valeurs propres' |
---|
[630] | 105 | |
---|
[1489] | 106 | DO ii = 1, 2 |
---|
| 107 | ! |
---|
| 108 | DO ij = 1, ip1jmp1 |
---|
| 109 | zu(ij) = RAN1(idum) -.5 |
---|
| 110 | ENDDO |
---|
| 111 | CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) |
---|
| 112 | DO ij = 1, ip1jm |
---|
| 113 | zv(ij) = RAN1(idum) -.5 |
---|
| 114 | ENDDO |
---|
| 115 | CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) |
---|
[630] | 116 | |
---|
[1489] | 117 | CALL minmax(iip1*jjp1,zu,umin,ullm ) |
---|
| 118 | CALL minmax(iip1*jjm, zv,vmin,vllm ) |
---|
[630] | 119 | |
---|
[1489] | 120 | ullm = ABS ( ullm ) |
---|
| 121 | vllm = ABS ( vllm ) |
---|
[630] | 122 | |
---|
[1489] | 123 | DO l = 1, 50 |
---|
| 124 | IF(ii.EQ.1) THEN |
---|
| 125 | !cccc CALL covcont( 1,zu,zv,zu,zv ) |
---|
| 126 | IF(lstardis) THEN |
---|
[1611] | 127 | CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy ) |
---|
[1489] | 128 | ELSE |
---|
[1611] | 129 | CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy ) |
---|
[1489] | 130 | ENDIF |
---|
| 131 | ELSE |
---|
| 132 | IF(lstardis) THEN |
---|
[1611] | 133 | CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy ) |
---|
[1489] | 134 | ELSE |
---|
[1611] | 135 | CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy ) |
---|
[1489] | 136 | ENDIF |
---|
| 137 | ENDIF |
---|
[630] | 138 | |
---|
[1611] | 139 | zllm = max(abs(maxval(gx)), abs(maxval(gy))) |
---|
| 140 | zu = gx / zllm |
---|
| 141 | zv = gy / zllm |
---|
[1489] | 142 | end DO |
---|
[630] | 143 | |
---|
[1489] | 144 | IF ( ii.EQ.1 ) THEN |
---|
| 145 | IF(lstardis) THEN |
---|
| 146 | cdivu = 1./zllm |
---|
| 147 | ELSE |
---|
| 148 | cdivu = zllm **( -1./nitergdiv ) |
---|
| 149 | ENDIF |
---|
| 150 | ELSE |
---|
| 151 | IF(lstardis) THEN |
---|
| 152 | crot = 1./ zllm |
---|
| 153 | ELSE |
---|
| 154 | crot = zllm **( -1./nitergrot ) |
---|
| 155 | ENDIF |
---|
| 156 | ENDIF |
---|
[630] | 157 | |
---|
[1489] | 158 | end DO |
---|
[630] | 159 | |
---|
[1489] | 160 | ! petit test pour les operateurs non star: |
---|
| 161 | ! ---------------------------------------- |
---|
[630] | 162 | |
---|
[1489] | 163 | ! IF(.NOT.lstardis) THEN |
---|
| 164 | fact = rad*24./REAL(jjm) |
---|
| 165 | fact = fact*fact |
---|
[1490] | 166 | write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu |
---|
| 167 | write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot |
---|
| 168 | write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh |
---|
[1489] | 169 | ! ENDIF |
---|
[630] | 170 | |
---|
[1489] | 171 | !----------------------------------------------------------------------- |
---|
| 172 | ! variation verticale du coefficient de dissipation: |
---|
| 173 | ! -------------------------------------------------- |
---|
[630] | 174 | |
---|
[1697] | 175 | if (vert_prof_dissip == 1) then |
---|
[1489] | 176 | do l=1,llm |
---|
| 177 | pseudoz=8.*log(preff/presnivs(l)) |
---|
| 178 | zvert(l)=1+ & |
---|
| 179 | (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. & |
---|
| 180 | *(dissip_factz-1.) |
---|
| 181 | enddo |
---|
| 182 | else |
---|
| 183 | DO l=1,llm |
---|
| 184 | zvert(l)=1. |
---|
| 185 | ENDDO |
---|
| 186 | fact=2. |
---|
| 187 | DO l = 1, llm |
---|
| 188 | zz = 1. - preff/presnivs(l) |
---|
| 189 | zvert(l)= fact -( fact-1.)/( 1.+zz*zz ) |
---|
| 190 | ENDDO |
---|
| 191 | endif |
---|
[630] | 192 | |
---|
| 193 | |
---|
[1490] | 194 | write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale' |
---|
[630] | 195 | |
---|
[1489] | 196 | tetamin = 1.e+6 |
---|
[630] | 197 | |
---|
[1489] | 198 | DO l=1,llm |
---|
| 199 | tetaudiv(l) = zvert(l)/tetagdiv |
---|
| 200 | tetaurot(l) = zvert(l)/tetagrot |
---|
| 201 | tetah(l) = zvert(l)/tetatemp |
---|
[630] | 202 | |
---|
[1489] | 203 | IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) |
---|
| 204 | IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) |
---|
| 205 | IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) |
---|
| 206 | ENDDO |
---|
[630] | 207 | |
---|
[1502] | 208 | ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def |
---|
| 209 | IF (dissip_period == 0) THEN |
---|
| 210 | dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod |
---|
| 211 | write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period |
---|
| 212 | dissip_period = MAX(iperiod,dissip_period) |
---|
| 213 | END IF |
---|
[630] | 214 | |
---|
[1502] | 215 | dtdiss = dissip_period * dtvr |
---|
| 216 | write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr |
---|
| 217 | |
---|
[1489] | 218 | DO l = 1,llm |
---|
[1490] | 219 | write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), & |
---|
[1489] | 220 | dtdiss*tetah(l) |
---|
| 221 | ENDDO |
---|
[630] | 222 | |
---|
[1489] | 223 | END SUBROUTINE inidissip |
---|