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