Changeset 1454 for LMDZ5/trunk/libf/dyn3dpar
- Timestamp:
- Nov 18, 2010, 1:01:24 PM (14 years ago)
- Location:
- LMDZ5/trunk
- Files:
-
- 18 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk
- Property svn:mergeinfo changed
/LMDZ5/branches/LMDZ5V1.0-dev (added) merged: 1436-1438,1441-1449,1452-1453
- Property svn:mergeinfo changed
-
LMDZ5/trunk/libf/dyn3dpar/academic.h
r774 r1454 1 1 ! 2 ! $ Header$2 ! $Id$ 3 3 ! 4 real tetarappel(ip1jmp1,llm),taurappel 5 common/academic/tetarappel,taurappel 4 common/academic/tetarappel,knewt_t,kfrict,knewt_g,clat4 5 real :: tetarappel(ip1jmp1,llm) 6 real :: knewt_t(llm) 7 real :: kfrict(llm) 8 real :: knewt_g 9 real :: clat4(ip1jmp1) -
LMDZ5/trunk/libf/dyn3dpar/addfi_p.F
r1146 r1454 1 1 ! 2 ! $ Header$2 ! $Id$ 3 3 ! 4 4 SUBROUTINE addfi_p(pdt, leapf, forward, … … 7 7 USE parallel 8 8 USE infotrac, ONLY : nqtot 9 USE control_mod, ONLY : planet_type 9 10 IMPLICIT NONE 10 11 c … … 154 155 c$OMP END MASTER 155 156 156 DO iq = 1, 2 157 if (planet_type=="earth") then 158 ! earth case, special treatment for first 2 tracers (water) 159 DO iq = 1, 2 157 160 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 158 161 DO k = 1,llm … … 163 166 ENDDO 164 167 c$OMP END DO NOWAIT 165 ENDDO166 167 DO iq = 3, nqtot168 ENDDO 169 170 DO iq = 3, nqtot 168 171 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 169 172 DO k = 1,llm … … 174 177 ENDDO 175 178 c$OMP END DO NOWAIT 176 ENDDO 179 ENDDO 180 else 181 ! general case, treat all tracers equally) 182 DO iq = 1, nqtot 183 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 184 DO k = 1,llm 185 DO j = ijb,ije 186 pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt 187 pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt ) 188 ENDDO 189 ENDDO 190 c$OMP END DO NOWAIT 191 ENDDO 192 endif ! of if (planet_type=="earth") 177 193 178 194 c$OMP MASTER -
LMDZ5/trunk/libf/dyn3dpar/advtrac_p.F
r1403 r1454 132 132 ccc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) 133 133 c 134 ENDIF 134 ENDIF ! of IF(iadvtr.EQ.0) 135 135 136 136 iadvtr = iadvtr+1 … … 266 266 cym ----> Revérifier lors de la parallélisation des autres schemas 267 267 268 cym call massbar_p(massem,massebx,masseby) 268 cym call massbar_p(massem,massebx,masseby) 269 269 270 270 call vlspltgen_p( q,iadv, 2., massem, wg , … … 452 452 c$OMP BARRIER 453 453 454 ijb=ij_begin 455 ije=ij_end 454 if (planet_type=="earth") then 455 456 ijb=ij_begin 457 ije=ij_end 456 458 457 459 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 458 DO l = 1, llm460 DO l = 1, llm 459 461 DO ij = ijb, ije 460 462 finmasse(ij,l) = p(ij,l) - p(ij,l+1) 461 463 ENDDO 462 ENDDO464 ENDDO 463 465 c$OMP END DO 464 466 465 CALL qminimum_p( q, 2, finmasse )467 CALL qminimum_p( q, 2, finmasse ) 466 468 467 469 c------------------------------------------------------------------ … … 496 498 c$OMP BARRIER 497 499 iadvtr=0 500 endif ! of if (planet_type=="earth") 498 501 ENDIF ! if iadvtr.EQ.iapp_tracvl 499 502 -
LMDZ5/trunk/libf/dyn3dpar/caladvtrac_p.F
r1403 r1454 8 8 * flxw, pk, iapptrac) 9 9 USE parallel 10 USE infotrac 11 USE control_mod 10 USE infotrac, ONLY : nqtot 11 USE control_mod, ONLY : iapp_tracvl,planet_type 12 12 c 13 13 IMPLICIT NONE … … 30 30 c ---------- 31 31 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm),masse(ip1jmp1,llm) 32 REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqtot),dq( ip1jmp1,llm,2 ) 32 REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqtot) 33 real :: dq( ip1jmp1,llm,nqtot) 33 34 REAL teta( ip1jmp1,llm),pk( ip1jmp1,llm) 34 35 REAL :: flxw(ip1jmp1,llm) -
LMDZ5/trunk/libf/dyn3dpar/comconst.h
r1279 r1454 5 5 ! INCLUDE comconst.h 6 6 7 COMMON/comconst/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl, & 8 & dtvr,daysec, & 7 COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl, & 8 & iflag_top_bound 9 COMMON/comconstr/dtvr,daysec, & 9 10 & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg & 10 11 & ,dissip_factz,dissip_deltaz,dissip_zref & 11 & ,iflag_top_bound,tau_top_bound 12 & ,tau_top_bound, & 13 & daylen,year_day,molmass 12 14 13 15 14 16 INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl 15 REAL dtvr,daysec 16 REAL pi,dtphys,dtdiss,rad,r,cpp,kappa 17 REAL cotot,unsim,g,omeg 17 REAL dtvr ! dynamical time step (in s) 18 REAL daysec !length (in s) of a standard day 19 REAL pi ! something like 3.14159.... 20 REAL dtphys ! (s) time step for the physics 21 REAL dtdiss ! (s) time step for the dissipation 22 REAL rad ! (m) radius of the planet 23 REAL r ! Gas constant R=8.31 J.K-1.mol-1 24 REAL cpp ! Cp 25 REAL kappa ! kappa=R/Cp 26 REAL cotot 27 REAL unsim ! = 1./iim 28 REAL g ! (m/s2) gravity 29 REAL omeg ! (rad/s) rotation rate of the planet 18 30 REAL dissip_factz,dissip_deltaz,dissip_zref 19 31 INTEGER iflag_top_bound 20 32 REAL tau_top_bound 33 REAL daylen ! length of solar day, in 'standard' day length 34 REAL year_day ! Number of standard days in a year 35 REAL molmass ! (g/mol) molar mass of the atmosphere 21 36 22 37 -
LMDZ5/trunk/libf/dyn3dpar/conf_gcm.F
r1403 r1454 578 578 offline = .FALSE. 579 579 CALL getin('offline',offline) 580 580 IF (offline .AND. adjust) THEN 581 WRITE(lunout,*) 582 & 'WARNING : option offline does not work with adjust=y :' 583 WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ', 584 & 'and fluxstokev.nc will not be created' 585 WRITE(lunout,*) 586 & 'only the file phystoke.nc will still be created ' 587 END IF 588 581 589 !Config Key = config_inca 582 590 !Config Desc = Choix de configuration de INCA … … 768 776 offline = .FALSE. 769 777 CALL getin('offline',offline) 778 IF (offline .AND. adjust) THEN 779 WRITE(lunout,*) 780 & 'WARNING : option offline does not work with adjust=y :' 781 WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ', 782 & 'and fluxstokev.nc will not be created' 783 WRITE(lunout,*) 784 & 'only the file phystoke.nc will still be created ' 785 END IF 770 786 771 787 !Config Key = config_inca -
LMDZ5/trunk/libf/dyn3dpar/ener.h
r774 r1454 1 1 ! 2 ! $ Header$2 ! $Id$ 3 3 ! 4 c----------------------------------------------------------------------- 5 c INCLUDE 'ener.h' 4 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre 5 ! veillez à n'utiliser que des ! pour les commentaires 6 ! et à bien positionner les & des lignes de continuation 7 ! (les placer en colonne 6 et en colonne 73) 8 ! 9 ! INCLUDE 'ener.h' 6 10 7 COMMON/ener/ang0,etot0,ptot0,ztot0,stot0, 8 * ang,etot,ptot,ztot,stot,rmsdpdt ,9 *rmsv,gtot(llmm1)11 COMMON/ener/ang0,etot0,ptot0,ztot0,stot0, & 12 & ang,etot,ptot,ztot,stot,rmsdpdt , & 13 & rmsv,gtot(llmm1) 10 14 11 REAL ang0,etot0,ptot0,ztot0,stot0, 12 sang,etot,ptot,ztot,stot,rmsdpdt,rmsv,gtot15 REAL ang0,etot0,ptot0,ztot0,stot0, & 16 & ang,etot,ptot,ztot,stot,rmsdpdt,rmsv,gtot 13 17 14 c-----------------------------------------------------------------------18 !----------------------------------------------------------------------- -
LMDZ5/trunk/libf/dyn3dpar/fluxstokenc_p.F
r1403 r1454 57 57 58 58 c AC initialisations 59 cympbarug(:,:) = 0.59 pbarug(:,:) = 0. 60 60 cym pbarvg(:,:,:) = 0. 61 61 cym wg(:,:) = 0. -
LMDZ5/trunk/libf/dyn3dpar/friction_p.F
r1403 r1454 6 6 USE parallel 7 7 USE control_mod 8 #ifdef CPP_IOIPSL 9 USE IOIPSL 10 #else 11 ! if not using IOIPSL, we still need to use (a local version of) getin 12 USE ioipsl_getincom 13 #endif 8 14 IMPLICIT NONE 9 15 10 c=======================================================================11 c 12 c 13 c Objet: 14 c ------ 15 c 16 c *********** 17 c Friction 18 c *********** 19 c 20 c=======================================================================16 !======================================================================= 17 ! 18 ! Friction for the Newtonian case: 19 ! -------------------------------- 20 ! 2 possibilities (depending on flag 'friction_type' 21 ! friction_type=0 : A friction that is only applied to the lowermost 22 ! atmospheric layer 23 ! friction_type=1 : Friction applied on all atmospheric layer (but 24 ! (default) with stronger magnitude near the surface; see 25 ! iniacademic.F) 26 !======================================================================= 21 27 22 28 #include "dimensions.h" … … 24 30 #include "comgeom2.h" 25 31 #include "comconst.h" 26 27 REAL pdt 32 #include "iniprint.h" 33 #include "academic.h" 34 35 ! arguments: 36 REAL,INTENT(out) :: ucov( iip1,jjp1,llm ) 37 REAL,INTENT(out) :: vcov( iip1,jjm,llm ) 38 REAL,INTENT(in) :: pdt ! time step 39 40 ! local variables: 28 41 REAL modv(iip1,jjp1),zco,zsi 29 42 REAL vpn,vps,upoln,upols,vpols,vpoln 30 43 REAL u2(iip1,jjp1),v2(iip1,jjm) 31 REAL ucov( iip1,jjp1,llm ),vcov( iip1,jjm,llm ) 32 INTEGER i,j 33 REAL cfric 34 parameter (cfric=1.e-5) 44 INTEGER i,j,l 45 REAL,PARAMETER :: cfric=1.e-5 46 LOGICAL,SAVE :: firstcall=.true. 47 INTEGER,SAVE :: friction_type=1 48 CHARACTER(len=20) :: modname="friction_p" 49 CHARACTER(len=80) :: abort_message 50 !$OMP THREADPRIVATE(firstcall,friction_type) 35 51 integer :: jjb,jje 36 52 37 53 !$OMP SINGLE 54 IF (firstcall) THEN 55 ! set friction type 56 call getin("friction_type",friction_type) 57 if ((friction_type.lt.0).or.(friction_type.gt.1)) then 58 abort_message="wrong friction type" 59 write(lunout,*)'Friction: wrong friction type',friction_type 60 call abort_gcm(modname,abort_message,42) 61 endif 62 firstcall=.false. 63 ENDIF 64 !$OMP END SINGLE COPYPRIVATE(friction_type,firstcall) 65 66 if (friction_type.eq.0) then ! friction on first layer only 67 !$OMP SINGLE 38 68 c calcul des composantes au carre du vent naturel 39 69 jjb=jj_begin … … 138 168 vcov(iip1,j,1)=vcov(1,j,1) 139 169 enddo 170 !$OMP END SINGLE 171 endif ! of if (friction_type.eq.0) 172 173 if (friction_type.eq.1) then 174 ! for ucov() 175 jjb=jj_begin 176 jje=jj_end 177 if (pole_nord) jjb=jj_begin+1 178 if (pole_sud) jje=jj_end-1 179 180 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 181 do l=1,llm 182 ucov(1:iip1,jjb:jje,l)=ucov(1:iip1,jjb:jje,l)* 183 & (1.-pdt*kfrict(l)) 184 enddo 185 !$OMP END DO NOWAIT 186 187 ! for vcoc() 188 jjb=jj_begin 189 jje=jj_end 190 if (pole_sud) jje=jj_end-1 191 192 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 193 do l=1,llm 194 vcov(1:iip1,jjb:jje,l)=vcov(1:iip1,jjb:jje,l)* 195 & (1.-pdt*kfrict(l)) 196 enddo 197 !$OMP END DO 198 endif ! of if (friction_type.eq.1) 140 199 141 200 RETURN -
LMDZ5/trunk/libf/dyn3dpar/gcm.F
r1403 r1454 276 276 endif 277 277 278 if (planet_type.eq."earth") then 279 #ifdef CPP_EARTH 278 ! if (planet_type.eq."earth") then 280 279 ! Load an Earth-format start file 281 280 CALL dynetat0("start.nc",vcov,ucov, 282 281 & teta,q,masse,ps,phis, time_0) 283 #else 284 ! SW model also has Earth-format start files 285 ! (but can be used without the CPP_EARTH directive) 286 if (iflag_phys.eq.0) then 287 CALL dynetat0("start.nc",vcov,ucov, 288 & teta,q,masse,ps,phis, time_0) 289 endif 290 #endif 291 endif ! of if (planet_type.eq."earth") 282 ! endif ! of if (planet_type.eq."earth") 283 292 284 c write(73,*) 'ucov',ucov 293 285 c write(74,*) 'vcov',vcov … … 494 486 #endif 495 487 496 if (planet_type.eq."earth") then 488 ! if (planet_type.eq."earth") then 489 ! Write an Earth-format restart file 497 490 CALL dynredem0_p("restart.nc", day_end, phis) 498 endif491 ! endif 499 492 500 493 ecripar = .TRUE. -
LMDZ5/trunk/libf/dyn3dpar/grid_noro.F
r1403 r1454 458 458 C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS 459 459 460 PARAMETER (ISMo=300,JSMo=200) 461 REAL X(IMAR,JMAR),XF(ISMo,JSMo) 460 REAL X(IMAR,JMAR),XF(IMAR,JMAR) 462 461 real WEIGHTpb(-1:1,-1:1) 463 462 464 if(imar.gt.ismo) stop'surdimensionner ismo dans mva9 (grid_noro)'465 if(jmar.gt.jsmo) stop'surdimensionner jsmo dans mva9 (grid_noro)'466 463 467 464 SUM=0. -
LMDZ5/trunk/libf/dyn3dpar/infotrac.F90
r1403 r1454 61 61 CHARACTER(len=1), DIMENSION(3) :: txts 62 62 CHARACTER(len=2), DIMENSION(9) :: txtp 63 CHARACTER(len= 13) :: str1,str263 CHARACTER(len=23) :: str1,str2 64 64 65 65 INTEGER :: nqtrue ! number of tracers read from tracer.def, without higer order of moment 66 66 INTEGER :: iq, new_iq, iiq, jq, ierr 67 INTEGER, EXTERNAL :: lnblnk 68 67 68 character(len=*),parameter :: modname="infotrac_init" 69 69 !----------------------------------------------------------------------- 70 70 ! Initialization : … … 100 100 OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr) 101 101 IF(ierr.EQ.0) THEN 102 WRITE(lunout,*) 'Open traceur.def : ok'102 WRITE(lunout,*) trim(modname),': Open traceur.def : ok' 103 103 READ(90,*) nqtrue 104 104 ELSE 105 WRITE(lunout,*) 'Problem in opening traceur.def' 106 WRITE(lunout,*) 'ATTENTION using defaut values' 107 nqtrue=4 ! Defaut value 108 END IF 109 ! Attention! Only for planet_type=='earth' 110 nbtr=nqtrue-2 105 WRITE(lunout,*) trim(modname),': Problem in opening traceur.def' 106 WRITE(lunout,*) trim(modname),': WARNING using defaut values' 107 if (planet_type=='earth') then 108 nqtrue=4 ! Default value for Earth 109 else 110 nqtrue=1 ! Default value for other planets 111 endif 112 END IF 113 if ( planet_type=='earth') then 114 ! For Earth, water vapour & liquid tracers are not in the physics 115 nbtr=nqtrue-2 116 else 117 ! Other planets (for now); we have the same number of tracers 118 ! in the dynamics than in the physics 119 nbtr=nqtrue 120 endif 111 121 ELSE 112 122 ! nbtr has been read from INCA by init_cont_lmdz() in gcm.F … … 114 124 END IF 115 125 116 IF ( nqtrue < 2) THEN117 WRITE(lunout,*) 'nqtrue=',nqtrue, ' is not allowded. 2 tracers is the minimum'126 IF ((planet_type=="earth").and.(nqtrue < 2)) THEN 127 WRITE(lunout,*) trim(modname),': nqtrue=',nqtrue, ' is not allowded. 2 tracers is the minimum' 118 128 CALL abort_gcm('infotrac_init','Not enough tracers',1) 119 129 END IF … … 156 166 ! Continue to read tracer.def 157 167 DO iq=1,nqtrue 158 READ(90, 999) hadv(iq),vadv(iq),tnom_0(iq)168 READ(90,*) hadv(iq),vadv(iq),tnom_0(iq) 159 169 END DO 160 170 CLOSE(90) 161 ELSE ! Without tracer.def 171 ELSE ! Without tracer.def, set default values 172 if (planet_type=="earth") then 173 ! for Earth, default is to have 4 tracers 162 174 hadv(1) = 14 163 175 vadv(1) = 14 … … 172 184 vadv(4) = 10 173 185 tnom_0(4) = 'PB' 186 else ! default for other planets 187 hadv(1) = 10 188 vadv(1) = 10 189 tnom_0(1) = 'dummy' 190 endif ! of if (planet_type=="earth") 174 191 END IF 175 192 176 WRITE(lunout,*) 'Valeur de traceur.def :'177 WRITE(lunout,*) 'nombre de traceurs ',nqtrue193 WRITE(lunout,*) trim(modname),': Valeur de traceur.def :' 194 WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue 178 195 DO iq=1,nqtrue 179 196 WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq) … … 217 234 new_iq=new_iq+10 ! 9 tracers added 218 235 ELSE 219 WRITE(lunout,*) 'This choice of advection schema is not available'236 WRITE(lunout,*) trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq) 220 237 CALL abort_gcm('infotrac_init','Bad choice of advection schema - 1',1) 221 238 END IF … … 227 244 nqtot = new_iq 228 245 229 WRITE(lunout,*) 'The choice of advection schema for one or more tracers'246 WRITE(lunout,*) trim(modname),': The choice of advection schema for one or more tracers' 230 247 WRITE(lunout,*) 'makes it necessary to add tracers' 231 WRITE(lunout,*) nqtrue,' is the number of true tracers'232 WRITE(lunout,*) nqtot, ' is the total number of tracers needed'248 WRITE(lunout,*) trim(modname)//': ',nqtrue,' is the number of true tracers' 249 WRITE(lunout,*) trim(modname)//': ',nqtot, ' is the total number of tracers needed' 233 250 234 251 ELSE … … 258 275 iadv(new_iq)=11 259 276 ELSE 260 WRITE(lunout,*)'This choice of advection schema is not available' 277 WRITE(lunout,*)trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq) 278 261 279 CALL abort_gcm('infotrac_init','Bad choice of advection schema - 2',1) 262 280 END IF … … 265 283 tname(new_iq)= tnom_0(iq) 266 284 IF (iadv(new_iq)==0) THEN 267 ttext(new_iq)= str1(1:lnblnk(str1))285 ttext(new_iq)=trim(str1) 268 286 ELSE 269 ttext(new_iq)= str1(1:lnblnk(str1))//descrq(iadv(new_iq))287 ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq)) 270 288 END IF 271 289 … … 276 294 new_iq=new_iq+1 277 295 iadv(new_iq)=-20 278 ttext(new_iq)= str2(1:lnblnk(str2))//txts(jq)279 tname(new_iq)= str1(1:lnblnk(str1))//txts(jq)296 ttext(new_iq)=trim(str2)//txts(jq) 297 tname(new_iq)=trim(str1)//txts(jq) 280 298 END DO 281 299 ELSE IF (iadv(new_iq)==30) THEN … … 283 301 new_iq=new_iq+1 284 302 iadv(new_iq)=-30 285 ttext(new_iq)= str2(1:lnblnk(str2))//txtp(jq)286 tname(new_iq)= str1(1:lnblnk(str1))//txtp(jq)303 ttext(new_iq)=trim(str2)//txtp(jq) 304 tname(new_iq)=trim(str1)//txtp(jq) 287 305 END DO 288 306 END IF … … 303 321 304 322 305 WRITE(lunout,*) 'Information stored in infotrac :'306 WRITE(lunout,*) 'iadv niadv tname ttext :'323 WRITE(lunout,*) trim(modname),': Information stored in infotrac :' 324 WRITE(lunout,*) trim(modname),': iadv niadv tname ttext :' 307 325 DO iq=1,nqtot 308 WRITE(lunout,*) iadv(iq),niadv(iq), tname(iq), ttext(iq) 326 WRITE(lunout,*) iadv(iq),niadv(iq),& 327 ' ',trim(tname(iq)),' ',trim(ttext(iq)) 309 328 END DO 310 329 … … 315 334 DO iq=1,nqtot 316 335 IF (iadv(iq)/=10 .AND. iadv(iq)/=14 .AND. iadv(iq)/=0) THEN 317 WRITE(lunout,*) 'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'336 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ' 318 337 CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1) 319 338 ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN 320 WRITE(lunout,*) 'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'339 WRITE(lunout,*)trim(modname),'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ' 321 340 CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1) 322 341 END IF … … 329 348 DEALLOCATE(tracnam) 330 349 331 999 FORMAT (i2,1x,i2,1x,a15)332 333 350 END SUBROUTINE infotrac_init 334 351 -
LMDZ5/trunk/libf/dyn3dpar/iniacademic.F
r1403 r1454 8 8 USE filtreg_mod 9 9 USE infotrac, ONLY : nqtot 10 USE control_mod 10 USE control_mod, ONLY: day_step,planet_type 11 #ifdef CPP_IOIPSL 12 USE IOIPSL 13 #else 14 ! if not using IOIPSL, we still need to use (a local version of) getin 15 USE ioipsl_getincom 16 #endif 17 USE Write_Field 11 18 12 19 … … 70 77 REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches 71 78 REAL phi(ip1jmp1,llm) ! geopotentiel 72 REAL ddsin,teta rappelj,tetarappell,zsig79 REAL ddsin,tetastrat,zsig,tetapv,w_pv ! variables auxiliaires 73 80 real tetajl(jjp1,llm) 74 81 INTEGER i,j,l,lsup,ij 75 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 76 89 real zz,ran1 77 90 integer idum … … 82 95 ! 1. Initializations for Earth-like case 83 96 ! -------------------------------------- 84 if (planet_type=="earth") then 85 c 97 c 98 ! initialize planet radius, rotation rate,... 99 call conf_planete 100 86 101 time_0=0. 87 day_ref= 0102 day_ref=1 88 103 annee_ref=0 89 104 90 105 im = iim 91 106 jm = jjm 92 day_ini = 0 93 omeg = 4.*asin(1.)/86400. 94 rad = 6371229. 95 g = 9.8 96 daysec = 86400. 107 day_ini = 1 97 108 dtvr = daysec/REAL(day_step) 98 109 zdtvr=dtvr 99 kappa = 0.2857143100 cpp = 1004.70885101 preff = 101325.102 pa = 50000.103 110 etot0 = 0. 104 111 ptot0 = 0. … … 120 127 if (.not.read_start) then 121 128 phis(:)=0. 122 q(:,:,1)=1.e-10 123 q(:,:,2)=1.e-15 124 q(:,:,3:nqtot)=0. 129 q(:,:,:)=0 125 130 CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) 126 131 endif … … 129 134 if (iflag_phys.eq.2) then 130 135 ! initializations for the academic case 131 ps(:)=1.e5 132 phis(:)=0. 133 c--------------------------------------------------------------------- 134 135 taurappel=10.*daysec 136 137 c--------------------------------------------------------------------- 138 c Calcul de la temperature potentielle : 139 c -------------------------------------- 140 136 137 ! if (planet_type=="earth") then 138 139 ! 1. local parameters 140 ! by convention, winter is in the southern hemisphere 141 ! Geostrophic wind or no wind? 142 ok_geost=.TRUE. 143 CALL getin('ok_geost',ok_geost) 144 ! Constants for Newtonian relaxation and friction 145 k_f=1. !friction 146 CALL getin('k_j',k_f) 147 k_f=1./(daysec*k_f) 148 k_c_s=4. !cooling surface 149 CALL getin('k_c_s',k_c_s) 150 k_c_s=1./(daysec*k_c_s) 151 k_c_a=40. !cooling free atm 152 CALL getin('k_c_a',k_c_a) 153 k_c_a=1./(daysec*k_c_a) 154 ! Constants for Teta equilibrium profile 155 teta0=315. ! mean Teta (S.H. 315K) 156 CALL getin('teta0',teta0) 157 ttp=200. ! Tropopause temperature (S.H. 200K) 158 CALL getin('ttp',ttp) 159 eps=0. ! Deviation to N-S symmetry(~0-20K) 160 CALL getin('eps',eps) 161 delt_y=60. ! Merid Temp. Gradient (S.H. 60K) 162 CALL getin('delt_y',delt_y) 163 delt_z=10. ! Vertical Gradient (S.H. 10K) 164 CALL getin('delt_z',delt_z) 165 ! Polar vortex 166 ok_pv=.false. 167 CALL getin('ok_pv',ok_pv) 168 phi_pv=-50. ! Latitude of edge of vortex 169 CALL getin('phi_pv',phi_pv) 170 phi_pv=phi_pv*pi/180. 171 dphi_pv=5. ! Width of the edge 172 CALL getin('dphi_pv',dphi_pv) 173 dphi_pv=dphi_pv*pi/180. 174 gam_pv=4. ! -dT/dz vortex (in K/km) 175 CALL getin('gam_pv',gam_pv) 176 177 ! 2. Initialize fields towards which to relax 178 ! Friction 179 knewt_g=k_c_a 141 180 DO l=1,llm 142 zsig=ap(l)/preff+bp(l) 143 if (zsig.gt.0.3) then 144 lsup=l 145 tetarappell=1./8.*(-log(zsig)-.5) 146 DO j=1,jjp1 147 ddsin=sin(rlatu(j))-sin(pi/20.) 148 tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell) 149 ENDDO 150 else 151 c Choix isotherme au-dessus de 300 mbar 152 do j=1,jjp1 153 tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa 154 enddo 155 endif ! of if (zsig.gt.0.3) 181 zsig=presnivs(l)/preff 182 knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) 183 kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) 184 ENDDO 185 DO j=1,jjp1 186 clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 187 ENDDO 188 189 ! Potential temperature 190 DO l=1,llm 191 zsig=presnivs(l)/preff 192 tetastrat=ttp*zsig**(-kappa) 193 tetapv=tetastrat 194 IF ((ok_pv).AND.(zsig.LT.0.1)) THEN 195 tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) 196 ENDIF 197 DO j=1,jjp1 198 ! Troposphere 199 ddsin=sin(rlatu(j)) 200 tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin 201 & -delt_z*(1.-ddsin*ddsin)*log(zsig) 202 ! Profil stratospherique isotherme (+vortex) 203 w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. 204 tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv 205 tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 206 ENDDO 156 207 ENDDO ! of DO l=1,llm 208 209 ! CALL writefield('theta_eq',tetajl) 157 210 158 211 do l=1,llm … … 165 218 enddo 166 219 167 c call dump2d(jjp1,llm,tetajl,'TEQ ') 168 169 CALL pression ( ip1jmp1, ap, bp, ps, p ) 170 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 171 CALL massdair(p,masse) 172 173 c intialisation du vent et de la temperature 174 teta(:,:)=tetarappel(:,:) 175 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) 176 call ugeostr(phi,ucov) 177 vcov=0. 178 q(:,:,1 )=1.e-10 179 q(:,:,2 )=1.e-15 180 q(:,:,3:nqtot)=0. 181 182 183 c perturbation aleatoire sur la temperature 184 idum = -1 185 zz = ran1(idum) 186 idum = 0 187 do l=1,llm 188 do ij=iip2,ip1jm 189 teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) 220 221 ! else 222 ! write(lunout,*)"iniacademic: planet types other than earth", 223 ! & " not implemented (yet)." 224 ! stop 225 ! endif ! of if (planet_type=="earth") 226 227 ! 3. Initialize fields (if necessary) 228 IF (.NOT. read_start) THEN 229 ! surface pressure 230 ps(:)=preff 231 ! ground geopotential 232 phis(:)=0. 233 234 CALL pression ( ip1jmp1, ap, bp, ps, p ) 235 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 236 CALL massdair(p,masse) 237 238 ! bulk initialization of temperature 239 teta(:,:)=tetarappel(:,:) 240 241 ! geopotential 242 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) 243 244 ! winds 245 if (ok_geost) then 246 call ugeostr(phi,ucov) 247 else 248 ucov(:,:)=0. 249 endif 250 vcov(:,:)=0. 251 252 ! bulk initialization of tracers 253 if (planet_type=="earth") then 254 ! Earth: first two tracers will be water 255 do i=1,nqtot 256 if (i.eq.1) q(:,:,i)=1.e-10 257 if (i.eq.2) q(:,:,i)=1.e-15 258 if (i.gt.2) q(:,:,i)=0. 259 enddo 260 else 261 q(:,:,:)=0 262 endif ! of if (planet_type=="earth") 263 264 ! add random perturbation to temperature 265 idum = -1 266 zz = ran1(idum) 267 idum = 0 268 do l=1,llm 269 do ij=iip2,ip1jm 270 teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) 271 enddo 190 272 enddo 191 enddo 192 193 do l=1,llm 194 do ij=1,ip1jmp1,iip1 195 teta(ij+iim,l)=teta(ij,l) 273 274 ! maintain periodicity in longitude 275 do l=1,llm 276 do ij=1,ip1jmp1,iip1 277 teta(ij+iim,l)=teta(ij,l) 278 enddo 196 279 enddo 197 enddo 198 199 200 201 c PRINT *,' Appel test_period avec tetarappel ' 202 c CALL test_period ( ucov,vcov,tetarappel,q,p,phis ) 203 c PRINT *,' Appel test_period avec teta ' 204 c CALL test_period ( ucov,vcov,teta,q,p,phis ) 205 206 c initialisation d'un traceur sur une colonne 207 j=jjp1*3/4 208 i=iip1/2 209 ij=(j-1)*iip1+i 210 q(ij,:,3)=1. 280 281 ENDIF ! of IF (.NOT. read_start) 211 282 endif ! of if (iflag_phys.eq.2) 212 283 213 else214 write(lunout,*)"iniacademic: planet types other than earth",215 & " not implemented (yet)."216 stop217 endif ! of if (planet_type=="earth")218 return219 284 END 220 285 c----------------------------------------------------------------------- -
LMDZ5/trunk/libf/dyn3dpar/initfluxsto_p.F
r1279 r1454 203 203 . llm, nivsigs, zvertiid) 204 204 c pour le fichier def 205 nivd(1) = 1 206 call histvert(filedid, 'sig_s', 'Niveaux sigma', 207 . 'sigma_level', 208 . 1, nivd, dvertiid) 209 205 if (mpi_rank==0) then 206 nivd(1) = 1 207 call histvert(filedid, 'sig_s', 'Niveaux sigma', 208 . 'sigma_level', 209 . 1, nivd, dvertiid) 210 endif 210 211 C 211 212 C Appels a histdef pour la definition des variables a sauvegarder … … 282 283 call histend(fileid) 283 284 call histend(filevid) 284 call histend(filedid)285 if (mpi_rank==0) call histend(filedid) 285 286 if (ok_sync) then 286 287 call histsync(fileid) 287 288 call histsync(filevid) 288 call histsync(filedid)289 if (mpi_rank==0) call histsync(filedid) 289 290 endif 290 291 -
LMDZ5/trunk/libf/dyn3dpar/integrd_p.F
r1403 r1454 6 6 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis,finvmaold) 7 7 USE parallel 8 USE control_mod 8 USE control_mod, only : planet_type 9 9 IMPLICIT NONE 10 10 … … 279 279 280 280 CALL qminimum_p( q, nq, deltap ) 281 endif ! of if (planet_type.eq."earth")282 281 c 283 282 c ..... Calcul de la valeur moyenne, unique aux poles pour q ..... … … 337 336 ENDDO 338 337 c$OMP END DO NOWAIT 338 339 endif ! of if (planet_type.eq."earth") 340 339 341 c 340 342 c -
LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F
r1403 r1454 234 234 235 235 c$OMP MASTER 236 dq =0.236 dq(:,:,:)=0. 237 237 CALL pression ( ip1jmp1, ap, bp, ps, p ) 238 238 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) … … 596 596 . flxw,pk, iapptrac) 597 597 598 IF (offline) THEN 599 Cmaf stokage du flux de masse pour traceurs OFF-LINE 600 601 #ifdef CPP_IOIPSL 602 CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis, 603 . dtvr, itau) 604 #endif 605 606 607 ENDIF ! of IF (offline) 608 c 598 C Stokage du flux de masse pour traceurs OFF-LINE 599 IF (offline .AND. .NOT. adjust) THEN 600 CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis, 601 . dtvr, itau) 602 ENDIF 603 609 604 ENDIF ! of IF( forward. OR . leapf ) 610 605 cc$OMP END PARALLEL … … 968 963 969 964 IF(iflag_phys.EQ.2) THEN ! "Newtonian" case 970 c Calcul academique de la physique = Rappel Newtonien + fritcion 971 c -------------------------------------------------------------- 972 cym teta(:,:)=teta(:,:) 973 cym s -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel 965 ! Academic case : Simple friction and Newtonan relaxation 966 ! ------------------------------------------------------- 974 967 ijb=ij_begin 975 968 ije=ij_end 976 969 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 977 970 do l=1,llm 978 teta(ijb:ije,l)=teta(ijb:ije,l) 979 & -iphysiq*dtvr*(teta(ijb:ije,l)-tetarappel(ijb:ije,l))/taurappel 980 enddo 971 teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr* 972 & (teta(ijb:ije,l)-tetarappel(ijb:ije,l))* 973 & (knewt_g+knewt_t(l)*clat4(ijb:ije)) 974 enddo ! of do l=1,llm 981 975 !$OMP END DO 982 976 … … 987 981 call WaitRequest(Request_Physic) 988 982 c$OMP BARRIER 989 !$OMP MASTER 990 call friction_p(ucov,vcov,iphysiq*dtvr) 991 !$OMP END MASTER 983 call friction_p(ucov,vcov,dtvr) 992 984 !$OMP BARRIER 985 986 ! Sponge layer (if any) 987 IF (ok_strato) THEN 988 ! set dufi,dvfi,... to zero 989 ijb=ij_begin 990 ije=ij_end 991 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 992 do l=1,llm 993 dufi(ijb:ije,l)=0 994 dtetafi(ijb:ije,l)=0 995 dqfi(ijb:ije,l,1:nqtot)=0 996 enddo 997 !$OMP END DO 998 !$OMP SINGLE 999 dpfi(ijb:ije)=0 1000 !$OMP END SINGLE 1001 ijb=ij_begin 1002 ije=ij_end 1003 if (pole_sud) ije=ije-iip1 1004 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 1005 do l=1,llm 1006 dvfi(ijb:ije,l)=0 1007 enddo 1008 !$OMP END DO 1009 1010 CALL top_bound_p(vcov,ucov,teta,masse,dufi,dvfi,dtetafi) 1011 CALL addfi_p( dtvr, leapf, forward , 1012 $ ucov, vcov, teta , q ,ps , 1013 $ dufi, dvfi, dtetafi , dqfi ,dpfi ) 1014 !$OMP BARRIER 1015 ENDIF ! of IF (ok_strato) 993 1016 ENDIF ! of IF(iflag_phys.EQ.2) 994 1017 … … 1465 1488 c$OMP MASTER 1466 1489 1467 if (planet_type.eq."earth") then1490 ! if (planet_type.eq."earth") then 1468 1491 ! Write an Earth-format restart file 1469 1492 CALL dynredem1_p("restart.nc",0.0, 1470 1493 & vcov,ucov,teta,q,masse,ps) 1471 endif ! of if (planet_type.eq."earth")1494 ! endif ! of if (planet_type.eq."earth") 1472 1495 1473 1496 ! CLOSE(99) … … 1658 1681 1659 1682 IF(itau.EQ.itaufin) THEN 1660 if (planet_type.eq."earth") then1683 ! if (planet_type.eq."earth") then 1661 1684 c$OMP MASTER 1662 1685 CALL dynredem1_p("restart.nc",0.0, 1663 1686 . vcov,ucov,teta,q,masse,ps) 1664 1687 c$OMP END MASTER 1665 endif ! of if (planet_type.eq."earth")1688 ! endif ! of if (planet_type.eq."earth") 1666 1689 ENDIF ! of IF(itau.EQ.itaufin) 1667 1690 -
LMDZ5/trunk/libf/dyn3dpar/limit_netcdf.F90
r1425 r1454 97 97 kappa = 0.2857143 98 98 cpp = 1004.70885 99 dtvr = daysec/ FLOAT(day_step)99 dtvr = daysec/REAL(day_step) 100 100 CALL inigeom 101 101 … … 265 265 266 266 DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug) 267 #endif268 ! of #ifdef CPP_EARTH269 267 270 268 … … 592 590 593 591 !--- Mid-months times 594 mid_months(1)=0.5* FLOAT(mnth(1))592 mid_months(1)=0.5*REAL(mnth(1)) 595 593 DO k=2,nm 596 mid_months(k)=mid_months(k-1)+0.5* FLOAT(mnth(k-1)+mnth(k))594 mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k)) 597 595 END DO 598 596 … … 626 624 !------------------------------------------------------------------------------- 627 625 626 #endif 627 ! of #ifdef CPP_EARTH 628 628 629 629 END SUBROUTINE limit_netcdf
Note: See TracChangeset
for help on using the changeset viewer.