Changeset 609 for codes/icosagcm/devel
- Timestamp:
- 11/21/17 15:09:25 (7 years ago)
- Location:
- codes/icosagcm/devel/src
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/diagnostics/check_conserve.f90
r584 r609 127 127 WRITE(*,'(A,6E12.4)') 'AAM_vel : time,old,new,dissip,dyn,phys', dt*it, AAM_vel_old, AAM_vel, & 128 128 AAM_vel_source(AAM_dissip), AAM_vel_source(AAM_dyn), AAM_vel_source(AAM_phys) 129 WRITE(*,'(A,6E12.4)') 'AAM_tot : time,old,new,dissip,dyn,phys', dt*it, & 130 AAM_vel_old+AAM_mass_old, AAM_vel+AAM_mass, & 131 AAM_vel_source(AAM_dissip)+AAM_mass_source(AAM_dissip), & 132 AAM_vel_source(AAM_dyn) +AAM_mass_source(AAM_dyn), & 133 AAM_vel_source(AAM_phys) +AAM_mass_source(AAM_phys) 129 134 AAM_mass_old=AAM_mass 130 135 AAM_vel_old=AAM_vel -
codes/icosagcm/devel/src/initial/etat0_heldsz.f90
r549 r609 245 245 DO i=ii_begin-1,ii_end+1 246 246 ij=(j-1)*iim+i 247 theta(ij,l)=theta(ij,l) - dt*(theta(ij,l)-theta_eq(ij,l))* &247 theta(ij,l)=theta(ij,l) - itau_physics*dt*(theta(ij,l)-theta_eq(ij,l))* & 248 248 (knewt_g+knewt_t(l)*COS(lat(ij))**4 ) 249 249 ENDDO … … 253 253 254 254 Do l=1,llm 255 u(:,l)=u(:,l)*(1.- dt*kfrict(l))255 u(:,l)=u(:,l)*(1.-itau_physics*dt*kfrict(l)) 256 256 END DO 257 257 -
codes/icosagcm/devel/src/time/timeloop_gcm.f90
r601 r609 375 375 per_step = elapsed/(it-itau0) 376 376 throughput = dt/per_step 377 total = per_step* (itaumax-itau0)377 total = per_step*itaumax 378 378 WRITE(*,'(A,I5,A,F6.2,A,I6)') 'Time spent (s):',INT(elapsed), & 379 379 ' -- ms/step : ', 1000*per_step, & -
codes/icosagcm/devel/src/vertical/disvert.f90
r531 r609 23 23 SUBROUTINE init_disvert 24 24 USE disvert_std_mod, ONLY: ap_std=>ap, bp_std=>bp, presnivs_std=>presnivs, init_disvert_std=>init_disvert 25 USE disvert_strato_mod, ONLY: ap_strato=>ap, bp_strato=>bp, presnivs_strato=>presnivs, init_disvert_strato=>init_disvert 25 USE disvert_strato_mod, ONLY: ap_strato=>ap, bp_strato=>bp, presnivs_strato=>presnivs, init_disvert_strato=>init_disvert, init_disvert_strato_custom 26 26 USE disvert_apbp_mod, ONLY: ap_apbp=>ap, bp_apbp=>bp, presnivs_apbp=>presnivs, init_disvert_apbp=>init_disvert 27 27 USE disvert_ncar_mod, ONLY: ap_ncar=>ap, bp_ncar=>bp, presnivs_ncar=>presnivs, init_disvert_ncar=>init_disvert … … 69 69 70 70 CALL init_disvert_strato 71 ap=>ap_strato 72 bp=>bp_strato 73 presnivs=>presnivs_strato 74 75 CASE('strato_custom') 76 77 CALL init_disvert_strato_custom 71 78 ap=>ap_strato 72 79 bp=>bp_strato -
codes/icosagcm/devel/src/vertical/disvert_strato.f90
r550 r609 1 1 MODULE disvert_strato_mod 2 2 USE icosa 3 USE mpipara 4 USE getin_mod 5 USE earth_const 6 IMPLICIT NONE 7 PRIVATE 8 3 9 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:) 4 10 !$OMP THREADPRIVATE(ap) … … 8 14 !$OMP THREADPRIVATE(presnivs) 9 15 16 PUBLIC :: ap,bp,presnivs, init_disvert, init_disvert_strato_custom 17 10 18 CONTAINS 11 19 12 20 SUBROUTINE init_disvert 13 USE icosa 14 USE getin_mod 15 IMPLICIT NONE 16 REAL(rstd) :: dsigmin 21 REAL(rstd) :: dsigmin, snorm, x, dsig(llm), sig(llm+1) 22 INTEGER :: l 23 17 24 ALLOCATE(ap(llm+1)) 18 25 ALLOCATE(bp(llm+1)) 19 26 ALLOCATE(presnivs(llm)) 20 27 21 28 dsigmin=1.0 ! Should be 0.3 for CMIP5 22 29 CALL getin('disvert_dsigmin', dsigmin) 23 CALL disvert(dsigmin,ap,bp,presnivs) 24 25 END SUBROUTINE init_disvert 26 27 28 SUBROUTINE disvert(dsigmin,ap,bp,presnivs) 29 USE icosa 30 USE mpipara 31 USE earth_const 32 IMPLICIT NONE 33 REAL(rstd) :: dsigmin 34 REAL(rstd),INTENT(OUT) :: ap(:) 35 REAL(rstd),INTENT(OUT) :: bp(:) 36 REAL(rstd),INTENT(OUT) :: presnivs(:) 37 38 REAL(rstd) :: dsig(llm) 39 REAL(rstd) :: sig(llm+1) 40 REAL(rstd) :: snorm 41 REAL(rstd) :: x 42 INTEGER :: l 43 30 44 31 snorm = 0. 45 32 DO l = 1, llm … … 84 71 ENDIF 85 72 86 END SUBROUTINE disvert 73 END SUBROUTINE init_disvert 74 75 SUBROUTINE init_disvert_strato_custom 76 REAL(rstd) :: vert_scale,vert_dzmin,vert_dzlow,vert_z0low,vert_dzmid,vert_z0mid,vert_h_mid,vert_dzhig,vert_z0hig,vert_h_hig 77 REAL(rstd) :: z, snorm, dsig(llm), sig(llm+1), sig0(llm+1), zz(llm+1) 78 INTEGER :: l 79 ALLOCATE(ap(llm+1)) 80 ALLOCATE(bp(llm+1)) 81 ALLOCATE(presnivs(llm)) 82 !=================================================================== 83 ! David Cugnet, François Lott, Lionel Guez, Ehouoarn Millour, Fredho 84 ! 2014/05 85 ! custumize strato distribution 86 ! Al the parameter are given in km assuming a given scalehigh 87 vert_scale=7. ! scale hight 88 vert_dzmin=0.02 ! width of first layer 89 vert_dzlow=1. ! dz in the low atmosphere 90 vert_z0low=8. ! height at which resolution recches dzlow 91 vert_dzmid=3. ! dz in the mid atmsophere 92 vert_z0mid=70. ! height at which resolution recches dzmid 93 vert_h_mid=20. ! width of the transition 94 vert_dzhig=11. ! dz in the high atmsophere 95 vert_z0hig=80. ! height at which resolution recches dz 96 vert_h_hig=20. ! width of the transition 97 !=================================================================== 98 99 call getin('vert_scale',vert_scale) 100 call getin('vert_dzmin',vert_dzmin) 101 call getin('vert_dzlow',vert_dzlow) 102 call getin('vert_z0low',vert_z0low) 103 CALL getin('vert_dzmid',vert_dzmid) 104 CALL getin('vert_z0mid',vert_z0mid) 105 call getin('vert_h_mid',vert_h_mid) 106 call getin('vert_dzhig',vert_dzhig) 107 call getin('vert_z0hig',vert_z0hig) 108 call getin('vert_h_hig',vert_h_hig) 109 110 !!!!! scaleheight=vert_scale ! for consistency with further computations 111 ! Build geometrical distribution 112 zz(1)=0. 113 DO l=1,llm 114 z=zz(l) 115 zz(l+1)=zz(l)+vert_dzmin+vert_dzlow*TANH(z/vert_z0low)+ & 116 (vert_dzmid-vert_dzlow)* & 117 (TANH((z-vert_z0mid)/vert_h_mid)-TANH((-vert_z0mid)/vert_h_mid)) & 118 +(vert_dzhig-vert_dzmid-vert_dzlow)* & 119 (TANH((z-vert_z0hig)/vert_h_hig)-TANH((-vert_z0hig)/vert_h_hig)) 120 ENDDO 121 122 123 !=================================================================== 124 ! Comment added Fredho 2014/05/18, Saint-Louis, Senegal 125 ! From approximate z to ap, bp, so that p=ap+bp*p0 and p/p0=exp(-z/H) 126 ! sig0 is p/p0 127 ! sig is an intermediate distribution introduce to estimate bp 128 ! 1. sig0=exp(-z/H) 129 ! 2. inversion of sig0=(1-pa/p0)*sig+(1-pa/p0)*exp(1-1/sig**2) 130 ! 3. bp=exp(1-1/sig**2) 131 ! 4. ap deduced from the combination of 2 and 3 so that sig0=ap/p0+bp 132 !=================================================================== 133 134 sig0(llm+1)=0. 135 DO l = 1, llm 136 sig0(l)=EXP(-zz(l)/vert_scale) 137 !!! sig(l)=racinesig(sig0(l)) 138 ENDDO 139 sig=racinesig(sig0) 140 141 bp(llm+1) = 0. 142 DO l = 2, llm 143 bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) 144 ap(l) = pa * ( sig(l) - bp(l) ) 145 ENDDO 146 bp(1)=1. 147 !! ap=pa*(sig-bp) 148 ap(1)=0. 149 ap(llm+1) = pa * ( sig(llm+1) - bp(llm+1) ) 150 DO l = 1, llm 151 presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) 152 ENDDO 153 154 ! tell the world about it 155 IF (is_mpi_root) THEN 156 !$OMP MASTER 157 WRITE(*,*) "ap()=",ap 158 WRITE(*,*) "bp()=",bp 159 WRITE(*,*) "Approximative mid-layer pressure, assuming a surface pressure preff=",preff," Pa" 160 WRITE(*,*) "and approximative mid-layer height, assuming an atmospheric scale height of ",scale_height/1000," (km)" 161 DO l=1,llm 162 WRITE(*,*) 'PRESNIVS(',l,')=',presnivs(l),' Z ~ ',log(preff/presnivs(l))*scale_height/1000, & 163 ' DZ ~ ',scale_height/1000*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10)) 164 ENDDO 165 !$OMP END MASTER 166 ENDIF 167 168 CONTAINS 169 170 FUNCTION racinesig(sig) RESULT(sg) 171 ! 172 !------------------------------------------------------------------------------- 173 ! Fredho 2014/05/18 174 ! Purpose: Search for s solving (Pa/Preff)*sg+(1-Pa/Preff)*EXP(1-1./sg**2)=s 175 ! Notes: Uses Newton Raphson search 176 !------------------------------------------------------------------------------- 177 ! Arguments: 178 REAL, INTENT(IN) :: sig(:) 179 REAL(rstd) :: sg(SIZE(sig)) 180 !------------------------------------------------------------------------------- 181 ! Local variables: 182 INTEGER :: it, ns, maxit 183 REAL(rstd) :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib 184 !------------------------------------------------------------------------------- 185 ns=SIZE(sig); maxit=100 186 c1=Pa/Preff; c2=1.-c1 187 DO l=2,ns-1 188 sg(l)=sig(l) 189 DO it=1,maxit 190 f1=exp(1-1./sg(l)**2)*(1.-c1) 191 sg(l)=sg(l)-(c1*sg(l)+f1-sig(l))/(c1+2*f1*sg(l)**(-3)) 192 ENDDO 193 ! print*,'SSSSIG ',sig(l),sg(l),c1*sg(l)+exp(1-1./sg(l)**2)*(1.-c1) 194 ENDDO 195 sg(1)=1.; sg(ns)=0. 196 197 END FUNCTION racinesig 198 199 END SUBROUTINE init_disvert_strato_custom 87 200 88 201 END MODULE disvert_strato_mod
Note: See TracChangeset
for help on using the changeset viewer.