Index: trunk/ICOSA_LMDZ/src/disvert_icosa_lmdz.f90 =================================================================== --- trunk/ICOSA_LMDZ/src/disvert_icosa_lmdz.f90 (revision 2452) +++ trunk/ICOSA_LMDZ/src/disvert_icosa_lmdz.f90 (revision 2452) @@ -0,0 +1,184 @@ +module disvert_icosa_lmdz_mod + +implicit none + +contains + + SUBROUTINE disvert_icosa_lmdz(ap,bp,presnivs) + ! Routine to compute ap(),bp() and presnivs() in the same way as done + ! in LMDZ for planets (see LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F) + + USE prec, ONLY : rstd + USE getin_mod, ONLY : getin + USE grid_param, ONLY: llm + USE earth_const, ONLY: pa,preff,scale_height + USE abort_mod, ONLY: dynamico_abort + USE mpipara, ONLY: is_mpi_root + USE omp_para, ONLY: omp_in_parallel + USE transfert_omp_mod, ONLY: bcast_omp + USE free_unit_mod, ONLY : free_unit + IMPLICIT NONE + + REAL(rstd),INTENT(OUT) :: ap(:) + REAL(rstd),INTENT(OUT) :: bp(:) + REAL(rstd),INTENT(OUT) :: presnivs(:) + + INTEGER :: unit + CHARACTER(len=255) :: filename + INTEGER :: l,ok + LOGICAL :: hybrid + REAL(rstd) :: scaleheight ! atmospheric scale height (km) from file z2sig.def + REAL(rstd) :: zsig(llm) + REAL(rstd) :: sig(llm+1) + REAL(rstd) :: newsig + + filename="z2sig.def" !file to read + + ! User may choose between hybrid (default) or sigma coordinates + hybrid=.true. ! default: use hybrid coordinates + CALL getin('hybrid',hybrid) + ! User must also provide "pa" if in hybrid case + if (hybrid) then + pa=-999 ! absurd initial value to trigger error if not set + CALL getin('pa',pa) + IF (pa<0) THEN + WRITE(*,*) "disvert_icosa_lmdz error: pa=",pa + WRITE(*,*) " set a correct value for the pa parameter" + CALL dynamico_abort("disvert_icosa_lmdz : wrong value for pa") + ENDIF + endif + +!$OMP MASTER + unit=free_unit() + OPEN(unit,file=filename,status="old",action="read",iostat=ok) + IF (ok/=0) THEN + WRITE(*,*) "disvert_icosa_lmdz error: input file ",trim(filename)," not found!" + CALL dynamico_abort("disvert_icosa_lmdz : could not open z2sig.def file" ) + ENDIF + ! read "scaleheight" on the first line + READ(unit,fmt=*,iostat=ok) scaleheight + IF (ok/=0) THEN + WRITE(*,*) "disvert_icosa_lmdz error: failed reading scaleheight" + CALL dynamico_abort("disvert_icosa_lmdz : could not read input file") + ENDIF + ! While at it check that this "scaleheight" matches the "scale_height" + ! from the def files + IF ((1.-(scaleheight*1000.)/scale_height)>1.e-6) THEN + WRITE(*,*) "disvert_icosa_lmdz error: significant mismatch between" + WRITE(*,*) " scale_height (m) from def file:",scale_height + WRITE(*,*) " and scaleheight (km) from z2sig.def:",scaleheight + CALL dynamico_abort("disvert_icosa_lmdz : atmospheric scale height issue") + ENDIF + ! read zsig() line by line, starting from the surface up to model top + DO l=1,llm + READ(unit,fmt=*,iostat=ok) zsig(l) + IF (ok/=0) THEN + WRITE(*,*) "disvert_icosa_lmdz error: failed reading zsig(l) for l=",l + CALL dynamico_abort("disvert_icosa_lmdz : could not read input file" ) + ENDIF + ENDDO + CLOSE(unit) + + ! build sig() + sig(1)=1. + DO l=2,llm + sig(l)=0.5*(exp(-zsig(l)/scaleheight)+exp(-zsig(l-1)/scaleheight)) + ENDDO + sig(llm+1)=0 ! here we asume 0 pressure at top of atmosphere + + ! build ap() and bp() using sig(), pa and preff + IF (hybrid) THEN ! use hybrid coordinates + DO l=1,llm + CALL sig_hybrid(sig(l),pa,preff,newsig) + bp(l)=exp(1.-1./(newsig**2)) + ap(l)=pa*(newsig-bp(l)) + ENDDO + ELSE ! use sigma coordinates + ap(1:llm)=0 + bp(1:llm)=sig(1:llm) + ENDIF + ap(llm+1)=0 + bp(llm+1)=0 +!$OMP END MASTER + + IF (omp_in_parallel()) THEN + CALL bcast_omp(ap) + CALL bcast_omp(bp) + ENDIF + + ! build presnivs(), approximative mid-layer pressures + DO l=1,llm + presnivs(l)=0.5*(ap(l)+bp(l)*preff+ap(l+1)+bp(l+1)*preff) + ENDDO + + ! tell the world about it + IF (is_mpi_root) THEN +!$OMP MASTER + WRITE(*,*) "Vertical discretization set by disvert_icosa_lmdz with hybrid=",hybrid + WRITE(*,*) "ap()=",ap + WRITE(*,*) "bp()=",bp + WRITE(*,*) "Approximative mid-layer pressure, assuming a surface pressure preff=",preff," Pa" + WRITE(*,*) "and approximative mid-layer height, assuming an atmospheric scale height of ",scale_height/1000," (km)" + DO l=1,llm + WRITE(*,*) 'PRESNIVS(',l,')=',presnivs(l),' Z ~ ',log(preff/presnivs(l))*scale_height/1000, & + ' DZ ~ ',scale_height/1000*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10)) + ENDDO +!$OMP END MASTER + ENDIF + + CONTAINS + + SUBROUTINE sig_hybrid(sig,pa,preff,newsig) + ! F. Forget 2002 + ! Iterative procedure to find "newsig" the sigma value for a layer + ! based on the target values provided in z2sig.def (sig) + ! but such that: + ! (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig + ! iterative solving of the equation uses the fact that: + ! where exp(1-1./x**2) become << x + ! x exp(1-1./x**2)/x + ! 1 1 + ! 0.68 0.5 + ! 0.5 1.E-1 + ! 0.391 1.E-2 + ! 0.333 1.E-3 + ! 0.295 1.E-4 + ! 0.269 1.E-5 + ! 0.248 1.E-6 + ! => thus on can assume newsig = sig*preff/pa if sig*preff/pa < 0.25 + + IMPLICIT NONE + REAL(rstd),INTENT(in) :: sig + REAL(rstd),INTENT(in) :: pa + REAL(rstd),INTENT(in) :: preff + REAL(rstd),INTENT(out) :: newsig + REAL(rstd) :: F,x1,x2 + INTEGER :: j + + newsig=sig + x1=0 + x2=1 + IF (sig.ge.1) THEN ! handle sigma=1 case + newsig=sig + ELSE IF (sig*preff/pa.ge.0.25) THEN + DO J=1,9999 ! Overkill max number of iterations + F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig + IF (F.gt.1) THEN + x2=newsig + newsig=(x1+newsig)*0.5 + ELSE + x1=newsig + newsig=(x2+newsig)*0.5 + ENDIF + ! Convergence test: stop if sig is approximated to within 0.01m + IF (abs(10.*log(F)).LT.1.E-5) EXIT + ENDDO + ELSE ! IF (sig*preff/pa.le.0.25) THEN + newsig=sig*preff/pa + ENDIF + END SUBROUTINE sig_hybrid + + END SUBROUTINE disvert_icosa_lmdz + + +end module disvert_icosa_lmdz_mod Index: trunk/ICOSA_LMDZ/src/icosa_lmdz.f90 =================================================================== --- trunk/ICOSA_LMDZ/src/icosa_lmdz.f90 (revision 2424) +++ trunk/ICOSA_LMDZ/src/icosa_lmdz.f90 (revision 2452) @@ -1,4 +1,9 @@ PROGRAM icosa_lmdz - USE icosa_init_mod + USE icosa_init_mod, ONLY: icosa_init + USE disvert_plugin_mod, ONLY: disvert_plugin + USE disvert_icosa_lmdz_mod, ONLY: disvert_icosa_lmdz + + ! set up plugins + disvert_plugin => disvert_icosa_lmdz CALL icosa_init