source: codes/icosagcm/trunk/src/vertical/disvert_std.f90

Last change on this file was 669, checked in by dubos, 6 years ago

trunk/disvert : set pa=preff/2 rather than hard-coded pa=50000
File size: 2.1 KB
Line 
1MODULE disvert_std_mod
2  USE icosa
3  IMPLICIT NONE
4  PRIVATE
5
6  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:)
7!$OMP THREADPRIVATE(ap)
8  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: bp(:)
9!$OMP THREADPRIVATE(bp)
10  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:)
11!$OMP THREADPRIVATE(presnivs)
12
13  PUBLIC :: init_disvert, ap, bp, presnivs
14
15CONTAINS
16
17  SUBROUTINE init_disvert
18    pa = .5*preff ! differs from 50000 if preff differs from 1e5
19    CALL getin('pa',pa)
20
21    ALLOCATE(ap(llm+1))
22    ALLOCATE(bp(llm+1))
23    ALLOCATE(presnivs(llm))
24
25    CALL disvert(ap,bp,presnivs)   
26
27  END SUBROUTINE init_disvert 
28
29
30  SUBROUTINE disvert(ap,bp,presnivs)
31  USE mpipara
32  REAL(rstd),INTENT(OUT) :: ap(:)
33  REAL(rstd),INTENT(OUT) :: bp(:)
34  REAL(rstd),INTENT(OUT) :: presnivs(:)
35 
36  REAL(rstd) :: dsig(llm)
37  REAL(rstd) :: sig(llm+1)
38  REAL(rstd) :: snorm
39  INTEGER :: l
40 
41    snorm  = 0.
42    DO l = 1, llm
43      dsig(l) = 1.0 + 7.0 * SIN( Pi*(l-0.5)/(llm+1) )**2
44      snorm = snorm + dsig(l)
45    ENDDO   
46   
47    DO l = 1, llm
48      dsig(l) = dsig(l)/snorm
49    ENDDO
50
51    sig(llm+1) = 0.
52    DO l = llm, 1, -1
53      sig(l) = sig(l+1) + dsig(l)
54    ENDDO
55
56    bp(llm+1) =   0.
57    DO l = 1, llm
58      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
59      ap(l) = pa * ( sig(l) - bp(l) )
60    ENDDO
61    bp(1)=1.
62    ap(1)=0.
63    ap(llm+1) = pa * ( sig(llm+1) - bp(llm+1) )
64    DO l = 1, llm
65      presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
66    ENDDO
67   
68    ! tell the world about it
69    IF (is_mpi_root) THEN
70!$OMP MASTER
71      WRITE(*,*) "ap()=",ap
72      WRITE(*,*) "bp()=",bp
73      WRITE(*,*) "Approximative mid-layer pressure, assuming a surface pressure preff=",preff," Pa"
74      WRITE(*,*) "and approximative mid-layer height, assuming an atmospheric scale height of ",scale_height/1000," (km)"
75      DO l=1,llm
76        WRITE(*,*) 'PRESNIVS(',l,')=',presnivs(l),'  Z ~ ',log(preff/presnivs(l))*scale_height/1000,       &
77                   ' DZ ~ ',scale_height/1000*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10))
78      ENDDO
79!$OMP END MASTER
80    ENDIF
81 
82  END SUBROUTINE disvert
83 
84END MODULE disvert_std_mod
Note: See TracBrowser for help on using the repository browser.