Changeset 66 for trunk/libf/dyn3d
- Timestamp:
- Feb 16, 2011, 4:57:45 PM (14 years ago)
- Location:
- trunk/libf/dyn3d
- Files:
-
- 3 edited
- 3 moved
Legend:
- Unmodified
- Added
- Removed
-
trunk/libf/dyn3d/caladvtrac.F
r7 r66 8 8 * flxw, pk) 9 9 c 10 USE infotrac 10 USE infotrac, ONLY : nqtot 11 11 USE control_mod, ONLY : iapp_tracvl,planet_type 12 12 -
trunk/libf/dyn3d/disvert.F90
r64 r66 1 ! 2 ! $Id: disvert.F 1403 2010-07-01 09:02:53Z fairhead $ 3 ! 4 SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) 1 ! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $ 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, *)'disvert: 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,*)'disvert: ATTENTION discretisation z a ajuster' 76 dsigmin=1. 77 endif 78 WRITE(LUNOUT,*) 'disvert: 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 = 2*asin(1.) * (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 *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 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, *)'disvert: BP ' 122 write(lunout, *) bp 123 write(lunout, *)'disvert: 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, *) 'disvert: 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 -
trunk/libf/dyn3d/etat0_netcdf.F90
r1 r66 1 1 ! 2 ! $Id: etat0_netcdf.F90 14 25 2010-09-02 13:45:23Z lguez$2 ! $Id: etat0_netcdf.F90 1486 2011-02-11 12:07:39Z fairhead $ 3 3 ! 4 4 !------------------------------------------------------------------------------- … … 98 98 REAL :: dummy 99 99 LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf 100 LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod 100 LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats 101 101 INTEGER :: iflag_radia, flag_aerosol 102 102 REAL :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut … … 130 130 !--- CONSTRUCT A GRID 131 131 CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & 132 callstats, & 132 133 solarlong0,seuil_inversion, & 133 134 fact_cldcon, facttemps,ok_newmicro,iflag_radia, & -
trunk/libf/dyn3d/infotrac.F90
r37 r66 177 177 ! Continue to read tracer.def 178 178 DO iq=1,nqtrue 179 READ(90, 999) hadv(iq),vadv(iq),tnom_0(iq)179 READ(90,*) hadv(iq),vadv(iq),tnom_0(iq) 180 180 END DO 181 181 CLOSE(90) … … 234 234 ! Continue to read tracer.def 235 235 DO iq=1,nqtrue 236 READ(90, 999) hadv(iq),vadv(iq),tnom_0(iq)236 READ(90,*) hadv(iq),vadv(iq),tnom_0(iq) 237 237 END DO 238 238 CLOSE(90) … … 316 316 tname(new_iq)= tnom_0(iq) 317 317 IF (iadv(new_iq)==0) THEN 318 ttext(new_iq)= str1(1:len_trim(str1))318 ttext(new_iq)=trim(str1) 319 319 ELSE 320 ttext(new_iq)= str1(1:len_trim(str1))//descrq(iadv(new_iq))320 ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq)) 321 321 END IF 322 322 … … 381 381 DEALLOCATE(tracnam) 382 382 383 999 FORMAT (i2,1x,i2,1x,a15)384 385 383 END SUBROUTINE infotrac_init 386 384 -
trunk/libf/dyn3d/iniacademic.F90
r64 r66 1 1 ! 2 ! $Id: iniacademic.F 1446 2010-10-22 09:27:25Z emillour$2 ! $Id: iniacademic.F90 1474 2011-01-14 11:04:45Z lguez $ 3 3 ! 4 c 5 c 6 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) 7 8 USE filtreg_mod 9 USE infotrac, ONLY : nqtot 10 USE control_mod, ONLY: day_step,planet_type 4 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) 5 6 USE filtreg_mod 7 USE infotrac, ONLY : nqtot 8 USE control_mod, ONLY: day_step,planet_type 11 9 #ifdef CPP_IOIPSL 12 10 USE IOIPSL 13 11 #else 14 ! if not using IOIPSL, we still need to use (a local version of) getin15 12 ! if not using IOIPSL, we still need to use (a local version of) getin 13 USE ioipsl_getincom 16 14 #endif 17 USE Write_Field 18 19 20 c%W% %G% 21 c======================================================================= 22 c 23 c Author: Frederic Hourdin original: 15/01/93 24 c ------- 25 c 26 c Subject: 27 c ------ 28 c 29 c Method: 30 c -------- 31 c 32 c Interface: 33 c ---------- 34 c 35 c Input: 36 c ------ 37 c 38 c Output: 39 c ------- 40 c 41 c======================================================================= 42 IMPLICIT NONE 43 c----------------------------------------------------------------------- 44 c Declararations: 45 c --------------- 46 47 #include "dimensions.h" 48 #include "paramet.h" 49 #include "comvert.h" 50 #include "comconst.h" 51 #include "comgeom.h" 52 #include "academic.h" 53 #include "ener.h" 54 #include "temps.h" 55 #include "iniprint.h" 56 #include "logic.h" 57 58 c Arguments: 59 c ---------- 60 61 real time_0 62 63 c variables dynamiques 64 REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants 65 REAL teta(ip1jmp1,llm) ! temperature potentielle 66 REAL q(ip1jmp1,llm,nqtot) ! champs advectes 67 REAL ps(ip1jmp1) ! pression au sol 68 REAL masse(ip1jmp1,llm) ! masse d'air 69 REAL phis(ip1jmp1) ! geopotentiel au sol 70 71 c Local: 72 c ------ 73 74 REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches 75 REAL pks(ip1jmp1) ! exner au sol 76 REAL pk(ip1jmp1,llm) ! exner au milieu des couches 77 REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches 78 REAL phi(ip1jmp1,llm) ! geopotentiel 79 REAL ddsin,tetastrat,zsig,tetapv,w_pv ! variables auxiliaires 80 real tetajl(jjp1,llm) 81 INTEGER i,j,l,lsup,ij 82 83 REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T 84 REAL k_f,k_c_a,k_c_s ! Constantes de rappel 85 LOGICAL ok_geost ! Initialisation vent geost. ou nul 86 LOGICAL ok_pv ! Polar Vortex 87 REAL phi_pv,dphi_pv,gam_pv ! Constantes pour polar vortex 88 89 real zz,ran1 90 integer idum 91 92 REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr 93 94 c----------------------------------------------------------------------- 95 ! 1. Initializations for Earth-like case 96 ! -------------------------------------- 97 c 98 99 print *, 'This is iniacademic' 100 101 ! initialize planet radius, rotation rate,... 102 call conf_planete 103 104 time_0=0. 105 day_ref=1 106 annee_ref=0 107 108 im = iim 109 jm = jjm 110 day_ini = 1 111 dtvr = daysec/REAL(day_step) 112 zdtvr=dtvr 113 etot0 = 0. 114 ptot0 = 0. 115 ztot0 = 0. 116 stot0 = 0. 117 ang0 = 0. 118 119 if (llm.eq.1) then 120 ! specific initializations for the shallow water case 121 kappa=1 15 USE Write_Field 16 17 ! Author: Frederic Hourdin original: 15/01/93 18 ! The forcing defined here is from Held and Suarez, 1994, Bulletin 19 ! of the American Meteorological Society, 75, 1825. 20 21 IMPLICIT NONE 22 23 ! Declararations: 24 ! --------------- 25 26 include "dimensions.h" 27 include "paramet.h" 28 include "comvert.h" 29 include "comconst.h" 30 include "comgeom.h" 31 include "academic.h" 32 include "ener.h" 33 include "temps.h" 34 include "iniprint.h" 35 include "logic.h" 36 37 ! Arguments: 38 ! ---------- 39 40 real time_0 41 42 ! variables dynamiques 43 REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants 44 REAL teta(ip1jmp1,llm) ! temperature potentielle 45 REAL q(ip1jmp1,llm,nqtot) ! champs advectes 46 REAL ps(ip1jmp1) ! pression au sol 47 REAL masse(ip1jmp1,llm) ! masse d'air 48 REAL phis(ip1jmp1) ! geopotentiel au sol 49 50 ! Local: 51 ! ------ 52 53 REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches 54 REAL pks(ip1jmp1) ! exner au sol 55 REAL pk(ip1jmp1,llm) ! exner au milieu des couches 56 REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches 57 REAL phi(ip1jmp1,llm) ! geopotentiel 58 REAL ddsin,zsig,tetapv,w_pv ! variables auxiliaires 59 real tetastrat ! potential temperature in the stratosphere, in K 60 real tetajl(jjp1,llm) 61 INTEGER i,j,l,lsup,ij 62 63 REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T 64 REAL k_f,k_c_a,k_c_s ! Constantes de rappel 65 LOGICAL ok_geost ! Initialisation vent geost. ou nul 66 LOGICAL ok_pv ! Polar Vortex 67 REAL phi_pv,dphi_pv,gam_pv ! Constantes pour polar vortex 68 69 real zz,ran1 70 integer idum 71 72 REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr 73 74 !----------------------------------------------------------------------- 75 ! 1. Initializations for Earth-like case 76 ! -------------------------------------- 77 ! 78 79 write(lunout,*) 'Iniacademic' 80 81 ! initialize planet radius, rotation rate,... 82 call conf_planete 83 84 time_0=0. 85 day_ref=1 86 annee_ref=0 87 88 im = iim 89 jm = jjm 90 day_ini = 1 91 dtvr = daysec/REAL(day_step) 92 zdtvr=dtvr 93 etot0 = 0. 94 ptot0 = 0. 95 ztot0 = 0. 96 stot0 = 0. 97 ang0 = 0. 98 99 if (llm == 1) then 100 ! specific initializations for the shallow water case 101 kappa=1 102 endif 103 104 CALL iniconst 105 CALL inigeom 106 CALL inifilr 107 108 if (llm == 1) then 109 ! initialize fields for the shallow water case, if required 110 if (.not.read_start) then 111 phis(:)=0. 112 q(:,:,:)=0 113 CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) 114 endif 115 endif 116 117 academic_case: if (iflag_phys == 2) then 118 ! initializations 119 120 ! 1. local parameters 121 ! by convention, winter is in the southern hemisphere 122 ! Geostrophic wind or no wind? 123 ok_geost=.TRUE. 124 CALL getin('ok_geost',ok_geost) 125 ! Constants for Newtonian relaxation and friction 126 k_f=1. !friction 127 CALL getin('k_j',k_f) 128 k_f=1./(daysec*k_f) 129 k_c_s=4. !cooling surface 130 CALL getin('k_c_s',k_c_s) 131 k_c_s=1./(daysec*k_c_s) 132 k_c_a=40. !cooling free atm 133 CALL getin('k_c_a',k_c_a) 134 k_c_a=1./(daysec*k_c_a) 135 ! Constants for Teta equilibrium profile 136 teta0=315. ! mean Teta (S.H. 315K) 137 CALL getin('teta0',teta0) 138 write(lunout,*) 'Iniacademic - teta0 ',teta0 139 write(lunout,*) 'Iniacademic - rad ',rad 140 ttp=200. ! Tropopause temperature (S.H. 200K) 141 CALL getin('ttp',ttp) 142 eps=0. ! Deviation to N-S symmetry(~0-20K) 143 CALL getin('eps',eps) 144 delt_y=60. ! Merid Temp. Gradient (S.H. 60K) 145 CALL getin('delt_y',delt_y) 146 delt_z=10. ! Vertical Gradient (S.H. 10K) 147 CALL getin('delt_z',delt_z) 148 ! Polar vortex 149 ok_pv=.false. 150 CALL getin('ok_pv',ok_pv) 151 phi_pv=-50. ! Latitude of edge of vortex 152 CALL getin('phi_pv',phi_pv) 153 phi_pv=phi_pv*pi/180. 154 dphi_pv=5. ! Width of the edge 155 CALL getin('dphi_pv',dphi_pv) 156 dphi_pv=dphi_pv*pi/180. 157 gam_pv=4. ! -dT/dz vortex (in K/km) 158 CALL getin('gam_pv',gam_pv) 159 160 ! 2. Initialize fields towards which to relax 161 ! Friction 162 knewt_g=k_c_a 163 DO l=1,llm 164 zsig=presnivs(l)/preff 165 knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) 166 kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) 167 ENDDO 168 DO j=1,jjp1 169 clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 170 ENDDO 171 172 ! Potential temperature 173 DO l=1,llm 174 zsig=presnivs(l)/preff 175 tetastrat=ttp*zsig**(-kappa) 176 tetapv=tetastrat 177 IF ((ok_pv).AND.(zsig.LT.0.1)) THEN 178 tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) 179 ENDIF 180 DO j=1,jjp1 181 ! Troposphere 182 ddsin=sin(rlatu(j)) 183 tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin & 184 -delt_z*(1.-ddsin*ddsin)*log(zsig) 185 !! Aymeric -- tests particuliers 186 if (planet_type=="giant") then 187 tetajl(j,l)=teta0+(delt_y* & 188 ((sin(rlatu(j)*3.14159*eps+0.0001))**2) & 189 / ((rlatu(j)*3.14159*eps+0.0001)**2)) & 190 -delt_z*log(zsig) 191 endif 192 ! Profil stratospherique isotherme (+vortex) 193 w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. 194 tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv 195 tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 196 ENDDO 197 ENDDO 198 199 ! CALL writefield('theta_eq',tetajl) 200 201 do l=1,llm 202 do j=1,jjp1 203 do i=1,iip1 204 ij=(j-1)*iip1+i 205 tetarappel(ij,l)=tetajl(j,l) 206 enddo 207 enddo 208 enddo 209 write(lunout,*) 'Iniacademic - check',tetajl(:,int(llm/2)),rlatu(:) 210 211 ! 3. Initialize fields (if necessary) 212 IF (.NOT. read_start) THEN 213 ! surface pressure 214 ps(:)=preff 215 ! ground geopotential 216 phis(:)=0. 217 218 CALL pression ( ip1jmp1, ap, bp, ps, p ) 219 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 220 CALL massdair(p,masse) 221 222 ! bulk initialization of temperature 223 teta(:,:)=tetarappel(:,:) 224 225 ! geopotential 226 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) 227 228 ! winds 229 if (ok_geost) then 230 call ugeostr(phi,ucov) 231 else 232 ucov(:,:)=0. 122 233 endif 123 124 CALL iniconst 125 CALL inigeom 126 CALL inifilr 127 128 if (llm.eq.1) then 129 ! initialize fields for the shallow water case, if required 130 if (.not.read_start) then 131 phis(:)=0. 132 q(:,:,:)=0 133 CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) 134 endif 135 endif 136 137 if (iflag_phys.eq.2) then 138 ! initializations for the academic case 139 140 ! if (planet_type=="earth") then 141 142 ! 1. local parameters 143 ! by convention, winter is in the southern hemisphere 144 ! Geostrophic wind or no wind? 145 ok_geost=.TRUE. 146 CALL getin('ok_geost',ok_geost) 147 ! Constants for Newtonian relaxation and friction 148 k_f=1. !friction 149 CALL getin('k_j',k_f) 150 k_f=1./(daysec*k_f) 151 k_c_s=4. !cooling surface 152 CALL getin('k_c_s',k_c_s) 153 k_c_s=1./(daysec*k_c_s) 154 k_c_a=40. !cooling free atm 155 CALL getin('k_c_a',k_c_a) 156 k_c_a=1./(daysec*k_c_a) 157 ! Constants for Teta equilibrium profile 158 teta0=315. ! mean Teta (S.H. 315K) 159 CALL getin('teta0',teta0) 160 print *, 'iniacademic - teta0 ', teta0 161 print *, 'iniacademic - rad ', rad 162 ttp=200. ! Tropopause temperature (S.H. 200K) 163 CALL getin('ttp',ttp) 164 eps=0. ! Deviation to N-S symmetry(~0-20K) 165 CALL getin('eps',eps) 166 delt_y=60. ! Merid Temp. Gradient (S.H. 60K) 167 CALL getin('delt_y',delt_y) 168 delt_z=10. ! Vertical Gradient (S.H. 10K) 169 CALL getin('delt_z',delt_z) 170 ! Polar vortex 171 ok_pv=.false. 172 CALL getin('ok_pv',ok_pv) 173 phi_pv=-50. ! Latitude of edge of vortex 174 CALL getin('phi_pv',phi_pv) 175 phi_pv=phi_pv*pi/180. 176 dphi_pv=5. ! Width of the edge 177 CALL getin('dphi_pv',dphi_pv) 178 dphi_pv=dphi_pv*pi/180. 179 gam_pv=4. ! -dT/dz vortex (in K/km) 180 CALL getin('gam_pv',gam_pv) 181 182 ! 2. Initialize fields towards which to relax 183 ! Friction 184 knewt_g=k_c_a 185 DO l=1,llm 186 zsig=presnivs(l)/preff 187 knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) 188 kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) 189 ENDDO 190 DO j=1,jjp1 191 clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 192 ENDDO 193 194 ! Potential temperature 195 DO l=1,llm 196 zsig=presnivs(l)/preff 197 tetastrat=ttp*zsig**(-kappa) 198 tetapv=tetastrat 199 IF ((ok_pv).AND.(zsig.LT.0.1)) THEN 200 tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) 201 ENDIF 202 DO j=1,jjp1 203 ! Troposphere 204 ddsin=sin(rlatu(j)) 205 tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin 206 & -delt_z*(1.-ddsin*ddsin)*log(zsig) 207 !! Aymeric -- tests particuliers 208 if (planet_type=="giant") then 209 tetajl(j,l)=teta0+(delt_y* 210 & ((sin(rlatu(j)*3.14159*eps+0.0001))**2) 211 & / ((rlatu(j)*3.14159*eps+0.0001)**2)) 212 & -delt_z*log(zsig) 213 !!! ddsin=sin(2.5*3.14159*rlatu(j)) 214 !!! tetajl(j,l)=teta0-delt_y*ddsin*ddsin 215 !!!! & -delt_z*(1.-ddsin*ddsin)*log(zsig) 216 endif 217 ! Profil stratospherique isotherme (+vortex) 218 w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. 219 tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv 220 tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 221 ENDDO 222 ENDDO ! of DO l=1,llm 223 224 ! CALL writefield('theta_eq',tetajl) 225 226 do l=1,llm 227 do j=1,jjp1 228 do i=1,iip1 229 ij=(j-1)*iip1+i 230 tetarappel(ij,l)=tetajl(j,l) 231 enddo 232 enddo 233 enddo 234 PRINT *, 'iniacademic - check',tetajl(:,int(llm/2)),rlatu(:) 235 236 237 ! else 238 ! write(lunout,*)"iniacademic: planet types other than earth", 239 ! & " not implemented (yet)." 240 ! stop 241 ! endif ! of if (planet_type=="earth") 242 243 ! 3. Initialize fields (if necessary) 244 IF (.NOT. read_start) THEN 245 ! surface pressure 246 ps(:)=preff 247 ! ground geopotential 248 phis(:)=0. 249 250 CALL pression ( ip1jmp1, ap, bp, ps, p ) 251 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 252 CALL massdair(p,masse) 253 254 ! bulk initialization of temperature 255 teta(:,:)=tetarappel(:,:) 256 257 ! geopotential 258 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) 259 260 ! winds 261 if (ok_geost) then 262 call ugeostr(phi,ucov) 263 else 264 ucov(:,:)=0. 265 endif 266 vcov(:,:)=0. 267 268 ! bulk initialization of tracers 269 if (planet_type=="earth") then 270 ! Earth: first two tracers will be water 271 do i=1,nqtot 272 if (i.eq.1) q(:,:,i)=1.e-10 273 if (i.eq.2) q(:,:,i)=1.e-15 274 if (i.gt.2) q(:,:,i)=0. 275 enddo 276 else 277 q(:,:,:)=0 278 endif ! of if (planet_type=="earth") 279 280 ! add random perturbation to temperature 281 idum = -1 282 zz = ran1(idum) 283 idum = 0 284 do l=1,llm 285 do ij=iip2,ip1jm 286 teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) 287 enddo 288 enddo 289 290 ! maintain periodicity in longitude 291 do l=1,llm 292 do ij=1,ip1jmp1,iip1 293 teta(ij+iim,l)=teta(ij,l) 294 enddo 295 enddo 296 297 ENDIF ! of IF (.NOT. read_start) 298 endif ! of if (iflag_phys.eq.2) 299 300 END 301 c----------------------------------------------------------------------- 234 vcov(:,:)=0. 235 236 ! bulk initialization of tracers 237 if (planet_type=="earth") then 238 ! Earth: first two tracers will be water 239 do i=1,nqtot 240 if (i == 1) q(:,:,i)=1.e-10 241 if (i == 2) q(:,:,i)=1.e-15 242 if (i.gt.2) q(:,:,i)=0. 243 enddo 244 else 245 q(:,:,:)=0 246 endif ! of if (planet_type=="earth") 247 248 ! add random perturbation to temperature 249 idum = -1 250 zz = ran1(idum) 251 idum = 0 252 do l=1,llm 253 do ij=iip2,ip1jm 254 teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) 255 enddo 256 enddo 257 258 ! maintain periodicity in longitude 259 do l=1,llm 260 do ij=1,ip1jmp1,iip1 261 teta(ij+iim,l)=teta(ij,l) 262 enddo 263 enddo 264 265 ENDIF ! of IF (.NOT. read_start) 266 endif academic_case 267 268 END SUBROUTINE iniacademic -
trunk/libf/dyn3d/ugeostr.F90
r64 r66 1 1 ! 2 ! $Id: ugeostr.F 1403 2010-07-01 09:02:53Z fairhead$2 ! $Id: ugeostr.F90 1474 2011-01-14 11:04:45Z lguez $ 3 3 ! 4 4 subroutine ugeostr(phi,ucov) 5 5 6 ! Calcul du vent covariant geostrophique a partir du champ de 7 ! geopotentiel. 8 ! We actually compute: (1 - cos^8 \phi) u_g 9 ! to have a wind going smoothly to 0 at the equator. 10 ! We assume that the surface pressure is uniform so that model 11 ! levels are pressure levels. 6 12 7 c Calcul du vent covariant geostrophique a partir du champs de 8 c geopotentiel en supposant que le vent au sol est nul. 13 implicit none 9 14 10 implicit none 15 include "dimensions.h" 16 include "paramet.h" 17 include "comconst.h" 18 include "comgeom2.h" 11 19 12 #include "dimensions.h" 13 #include "paramet.h" 14 #include "comconst.h" 15 #include "comgeom2.h" 20 real ucov(iip1,jjp1,llm),phi(iip1,jjp1,llm) 21 real um(jjm,llm),fact,u(iip1,jjm,llm) 22 integer i,j,l 16 23 17 real ucov(iip1,jjp1,llm),phi(iip1,jjp1,llm) 18 real um(jjm,llm),fact,u(iip1,jjm,llm) 19 integer i,j,l 24 real zlat 20 25 21 real zlat26 um(:,:)=0 ! initialize um() 22 27 23 um(:,:)=0 ! initialize um()28 DO j=1,jjm 24 29 25 DO j=1,jjm 30 if (abs(sin(rlatv(j))).lt.1.e-4) then 31 zlat=1.e-4 32 else 33 zlat=rlatv(j) 34 endif 35 fact=cos(zlat) 36 fact=fact*fact 37 fact=fact*fact 38 fact=fact*fact 39 fact=(1.-fact)/ & 40 (2.*omeg*sin(zlat)*(rlatu(j+1)-rlatu(j))) 41 fact=-fact/rad 42 DO l=1,llm 43 DO i=1,iim 44 u(i,j,l)=fact*(phi(i,j+1,l)-phi(i,j,l)) 45 um(j,l)=um(j,l)+u(i,j,l)/REAL(iim) 46 ENDDO 47 ENDDO 48 ENDDO 49 call dump2d(jjm,llm,um,'Vent-u geostrophique') 26 50 27 if (abs(sin(rlatv(j))).lt.1.e-4) then 28 zlat=1.e-4 29 else 30 zlat=rlatv(j) 31 endif 32 fact=cos(zlat) 33 fact=fact*fact 34 fact=fact*fact 35 fact=fact*fact 36 fact=(1.-fact)/ 37 s (2.*omeg*sin(zlat)*(rlatu(j+1)-rlatu(j))) 38 fact=-fact/rad 39 DO l=1,llm 40 DO i=1,iim 41 u(i,j,l)=fact*(phi(i,j+1,l)-phi(i,j,l)) 42 um(j,l)=um(j,l)+u(i,j,l)/REAL(iim) 43 ENDDO 44 ENDDO 45 ENDDO 46 call dump2d(jjm,llm,um,'Vent-u geostrophique') 51 ! calcul des champ de vent: 47 52 48 c 49 c----------------------------------------------------------------------- 50 c calcul des champ de vent: 51 c ------------------------- 53 DO l=1,llm 54 DO i=1,iip1 55 ucov(i,1,l)=0. 56 ucov(i,jjp1,l)=0. 57 end DO 58 DO j=2,jjm 59 DO i=1,iim 60 ucov(i,j,l) = 0.5*(u(i,j,l)+u(i,j-1,l))*cu(i,j) 61 end DO 62 ucov(iip1,j,l)=ucov(1,j,l) 63 end DO 64 end DO 52 65 53 DO 301 l=1,llm 54 DO 302 i=1,iip1 55 ucov(i,1,l)=0. 56 ucov(i,jjp1,l)=0. 57 302 CONTINUE 58 DO 304 j=2,jjm 59 DO 305 i=1,iim 60 ucov(i,j,l) = 0.5*(u(i,j,l)+u(i,j-1,l))*cu(i,j) 61 305 CONTINUE 62 ucov(iip1,j,l)=ucov(1,j,l) 63 304 CONTINUE 64 301 CONTINUE 66 print *, 301 65 67 66 print*,301 67 68 return 69 end 68 end subroutine ugeostr
Note: See TracChangeset
for help on using the changeset viewer.