source: codes/icosagcm/trunk/src/vertical/disvert_strato.f90 @ 606

Last change on this file since 606 was 606, checked in by dubos, 7 years ago

trunk : new disvert strato_custom
File size: 6.9 KB
Line 
1MODULE disvert_strato_mod
2  USE icosa
3  USE mpipara
4  USE getin_mod
5  USE earth_const
6  IMPLICIT NONE
7  PRIVATE
8
9  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:)
10!$OMP THREADPRIVATE(ap)
11  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: bp(:)
12!$OMP THREADPRIVATE(bp)
13  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:)
14!$OMP THREADPRIVATE(presnivs)
15
16  PUBLIC :: ap,bp,presnivs, init_disvert, init_disvert_strato_custom
17
18CONTAINS
19 
20  SUBROUTINE init_disvert
21    REAL(rstd) :: dsigmin, snorm, x, dsig(llm),  sig(llm+1)
22    INTEGER :: l
23
24    ALLOCATE(ap(llm+1))
25    ALLOCATE(bp(llm+1))
26    ALLOCATE(presnivs(llm))
27   
28    dsigmin=1.0  ! Should be 0.3 for CMIP5
29    CALL getin('disvert_dsigmin', dsigmin) 
30   
31    snorm  = 0.
32    DO l = 1, llm
33      x = Pi*(l-0.5)/(llm+1)
34      dsig(l) = (dsigmin + 7.0 * SIN(x)**2) * (0.5* (1.-tanh( (x-Pi/2)/(Pi/2) ) ))**2
35      snorm = snorm + dsig(l)
36    ENDDO   
37   
38    DO l = 1, llm
39      dsig(l) = dsig(l)/snorm
40    ENDDO
41
42    sig(llm+1) = 0.
43    DO l = llm, 1, -1
44      sig(l) = sig(l+1) + dsig(l)
45    ENDDO
46
47    bp(llm+1) =   0.
48    DO l = 1, llm
49      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
50      ap(l) = pa * ( sig(l) - bp(l) )
51    ENDDO
52    bp(1)=1.
53    ap(1)=0.
54    ap(llm+1) = pa * ( sig(llm+1) - bp(llm+1) )
55    DO l = 1, llm
56      presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
57    ENDDO
58   
59    ! tell the world about it
60    IF (is_mpi_root) THEN
61!$OMP MASTER
62      WRITE(*,*) "ap()=",ap
63      WRITE(*,*) "bp()=",bp
64      WRITE(*,*) "Approximative mid-layer pressure, assuming a surface pressure preff=",preff," Pa"
65      WRITE(*,*) "and approximative mid-layer height, assuming an atmospheric scale height of ",scale_height/1000," (km)"
66      DO l=1,llm
67        WRITE(*,*) 'PRESNIVS(',l,')=',presnivs(l),'  Z ~ ',log(preff/presnivs(l))*scale_height/1000,       &
68                   ' DZ ~ ',scale_height/1000*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10))
69      ENDDO
70!$OMP END MASTER
71    ENDIF
72 
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
200 
201END MODULE disvert_strato_mod
Note: See TracBrowser for help on using the repository browser.