source: codes/icosagcm/trunk/src/disvert.f90 @ 176

Last change on this file since 176 was 166, checked in by dubos, 11 years ago

New flag disvert=none
File size: 7.3 KB
Line 
1MODULE disvert_mod
2  USE icosa
3  REAL(rstd), SAVE, POINTER :: ap(:)
4  REAL(rstd), SAVE, POINTER :: bp(:)
5  REAL(rstd), SAVE, POINTER :: presnivs(:)
6  REAL(rstd), SAVE, POINTER :: mass_al(:), mass_bl(:), mass_ak(:), mass_bk(:), mass_dak(:), mass_dbk(:)
7
8  REAL(rstd) :: ptop ! pressure at top of atmosphere l=llm+1
9  ! vertical coordinate : Lagrangian or mass-based
10  INTEGER, SAVE :: caldyn_eta
11  INTEGER, PARAMETER :: eta_mass=1, eta_lag=2
12  LOGICAL,SAVE :: ap_bp_present
13
14CONTAINS
15
16  SUBROUTINE init_disvert
17  USE disvert_std_mod, ONLY: ap_std=>ap, bp_std=>bp, presnivs_std=>presnivs, init_disvert_std=>init_disvert
18  USE disvert_ncar_mod, ONLY: ap_ncar=>ap, bp_ncar=>bp, presnivs_ncar=>presnivs, init_disvert_ncar=>init_disvert
19  USE disvert_ncarl30_mod, ONLY: ap_ncarl30=>ap, bp_ncarl30=>bp, presnivs_ncarl30=>presnivs, init_disvert_ncarl30=>init_disvert
20  USE disvert_dcmip31_mod, ONLY: ap_dcmip31=>ap, bp_dcmip31=>bp, presnivs_dcmip31=>presnivs, init_disvert_dcmip31=>init_disvert
21  USE disvert_dcmip200_mod, ONLY: ap_dcmip200=>ap, bp_dcmip200=>bp, presnivs_dcmip200=>presnivs, init_disvert_dcmip200=>init_disvert
22  USE icosa
23  USE mpipara
24  IMPLICIT NONE
25    CHARACTER(len=255) :: def
26    INTEGER :: l
27
28    def='eta_mass'
29    CALL getin('caldyn_eta',def)
30    SELECT CASE(TRIM(def))
31    CASE('eta_mass')
32       caldyn_eta=eta_mass
33    CASE('eta_lag')
34       caldyn_eta=eta_lag
35    CASE DEFAULT
36       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_eta : <', TRIM(def),'> options are <eta_mass>, <eta_lag>'
37       STOP
38    END SELECT
39    IF (is_mpi_root) PRINT *, 'caldyn_eta=',TRIM(def), caldyn_eta, eta_lag, eta_mass
40
41    def='std'
42    ap_bp_present=.TRUE.
43    CALL getin("disvert",def)
44    SELECT CASE (TRIM(def))
45      CASE('none')
46         IF(caldyn_eta/=eta_lag) THEN
47            IF(is_mpi_root) PRINT *,'disvert_type=none must be used with caldyn_eta=eta_lag'
48            STOP
49         END IF
50         ap_bp_present=.FALSE.
51      CASE('std')
52   
53        CALL init_disvert_std
54        ap=>ap_std
55        bp=>bp_std
56        presnivs=>presnivs_std
57     
58      CASE ('ncar')
59
60        CALL init_disvert_ncar
61        ap=>ap_ncar
62        bp=>bp_ncar
63        presnivs=>presnivs_ncar
64
65      CASE ('dcmip31')
66
67        CALL init_disvert_dcmip31
68        ap=>ap_dcmip31
69        bp=>bp_dcmip31
70        presnivs=>presnivs_dcmip31
71
72      CASE ('dcmip200')
73
74        CALL init_disvert_dcmip200
75        ap=>ap_dcmip200
76        bp=>bp_dcmip200
77        presnivs=>presnivs_dcmip200
78
79      CASE ('ncarl30')
80
81        CALL init_disvert_ncarl30
82        ap=>ap_ncarl30
83        bp=>bp_ncarl30
84        presnivs=>presnivs_ncarl30
85       
86      CASE default
87        IF (is_mpi_root) PRINT*,'Bad selector for variable disvert : <', TRIM(def),"> options are <std>, <ncar>, <ncarl30>" 
88        STOP
89       
90    END SELECT
91
92    IF(ap_bp_present) THEN
93       ! Convert from pressure coordinate to mass coordinate
94       ! pk = ap(k) + bp(k)*ps(ij) = ptop + g*( mass_ak(k) + mass_bk(k)*ms(ij) )
95       ptop = ap(llm+1)
96       ALLOCATE(mass_al(llm+1))
97       ALLOCATE(mass_bl(llm+1))
98       ALLOCATE(mass_ak(llm))
99       ALLOCATE(mass_bk(llm))
100       ALLOCATE(mass_dak(llm))
101       ALLOCATE(mass_dbk(llm))
102
103       ! FIXME : leave ps for the moment ; change ps to ms later
104       DO l=1,llm+1
105          !       mass_al(l) = (ap(l)-ptop)/g
106          !       mass_bl(l) = bp(l)/g
107          mass_al(l) = (ap(l)-ptop)
108          mass_bl(l) = bp(l)
109       END DO
110       DO l=1,llm
111          mass_ak(l) = .5*(mass_al(l)+mass_al(l+1))
112          mass_bk(l) = .5*(mass_bl(l)+mass_bl(l+1))
113          mass_dak(l) = mass_al(l)-mass_al(l+1)  ! >0
114          mass_dbk(l) = mass_bl(l)-mass_bl(l+1)  ! >0
115       END DO
116    END IF
117
118  END SUBROUTINE init_disvert 
119 
120  SUBROUTINE compute_rhodz(comp, ps, rhodz)
121    USE icosa
122    USE omp_para
123    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
124    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
125    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
126    REAL(rstd) :: m, err
127    INTEGER :: l,i,j,ij,dd
128    err=0.
129
130    IF(comp) THEN
131       dd=1
132    ELSE
133       dd=0
134    END IF
135
136    DO l = ll_begin, ll_end
137       DO j=jj_begin-dd,jj_end+dd
138          DO i=ii_begin-dd,ii_end+dd
139             ij=(j-1)*iim+i
140             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
141             IF(comp) THEN
142                rhodz(ij,l) = m
143             ELSE
144                err = MAX(err,abs(m-rhodz(ij,l)))
145             END IF
146          ENDDO
147       ENDDO
148    ENDDO
149
150    IF(.NOT. comp) THEN
151       IF(err>1e-10) THEN
152          PRINT *, 'Discrepancy between ps and rhodz detected', err
153          STOP
154       ELSE
155!          PRINT *, 'No discrepancy between ps and rhodz detected'
156       END IF
157    END IF
158
159  END SUBROUTINE compute_rhodz
160
161 
162  SUBROUTINE write_apbp
163  USE icosa
164  USE netcdf_mod
165  IMPLICIT NONE
166    REAL(rstd) :: val(llm)
167    INTEGER :: status
168    INTEGER :: lev,ilev
169    INTEGER :: ncid,levid,ilevid,hyaiid,hybiid,hyamid,hybmid,P0id
170    INTEGER :: l
171   
172    IF(ap_bp_present) THEN
173    status = NF90_CREATE('apbp.nc', NF90_CLOBBER, ncid)
174    status = NF90_DEF_DIM(ncid,'lev',llm,lev)
175    status = NF90_DEF_DIM(ncid,'ilev',llm+1,ilev)
176   
177    status = NF90_DEF_VAR(ncid,'lev',NF90_DOUBLE,(/ lev /),levid)
178    status = NF90_PUT_ATT(ncid,levid,"long_name","hybrid level at midpoints (1000*(A+B))")
179    status = NF90_PUT_ATT(ncid,levid,"units","Pa")
180    status = NF90_PUT_ATT(ncid,levid,"positive","up")
181    status = NF90_PUT_ATT(ncid,levid,"standard_name","atmosphere_hybrid_sigma_pressure_coordinate")
182    status = NF90_PUT_ATT(ncid,levid,"formula_terms","a: hyam b: hybm p0: P0 ps: PS")
183
184    status = NF90_DEF_VAR(ncid,'ilev',NF90_DOUBLE,(/ ilev /),ilevid)
185    status = NF90_PUT_ATT(ncid,ilevid,"long_name","hybrid level at interfaces (1000*(A+B))")
186    status = NF90_PUT_ATT(ncid,ilevid,"units","Pa")
187    status = NF90_PUT_ATT(ncid,ilevid,"positive","up")
188    status = NF90_PUT_ATT(ncid,ilevid,"standard_name","atmosphere_hybrid_sigma_pressure_coordinate")
189    status = NF90_PUT_ATT(ncid,ilevid,"formula_terms","a: hyai b: hybi p0: P0 ps: PS")
190   
191    status = NF90_DEF_VAR(ncid,'hyai',NF90_DOUBLE,(/ ilev /),hyaiid)
192    status = NF90_PUT_ATT(ncid,hyaiid,"long_name","hybrid A coefficient at layer interfaces")
193
194    status = NF90_DEF_VAR(ncid,'hybi',NF90_DOUBLE,(/ ilev /),hybiid)
195    status = NF90_PUT_ATT(ncid,hybiid,"long_name","hybrid B coefficient at layer interfaces")
196
197    status = NF90_DEF_VAR(ncid,'hyam',NF90_DOUBLE,(/ lev /),hyamid)
198    status = NF90_PUT_ATT(ncid,hyamid,"long_name","hybrid A coefficient at midpoint interfaces")
199
200    status = NF90_DEF_VAR(ncid,'hybm',NF90_DOUBLE,(/ lev /),hybmid)
201    status = NF90_PUT_ATT(ncid,hybmid,"long_name","hybrid B coefficient at midpoint interfaces")
202   
203    status = NF90_DEF_VAR(ncid,'P0',NF90_DOUBLE,varid=P0id)
204
205    status = NF90_ENDDEF(ncid)   
206   
207    status=NF90_PUT_VAR(ncid,ilevid, ap(:)+bp(:)*Preff)
208   
209    DO l=1,llm
210      val(l)= 0.5*(ap(l+1)+ap(l)+Preff*(bp(l)+bp(l+1)))
211    ENDDO
212   
213    status=NF90_PUT_VAR(ncid,levid, val)
214
215    status=NF90_PUT_VAR(ncid,hyaiid, ap(:)/Preff)
216    status=NF90_PUT_VAR(ncid,hybiid, bp(:))
217   
218    DO l=1,llm
219      val(l)= 0.5*(ap(l+1)+ap(l))
220    ENDDO
221    status=NF90_PUT_VAR(ncid,hyamid, val(:)/Preff)
222   
223     DO l=1,llm
224      val(l)= 0.5*(bp(l+1)+bp(l))
225    ENDDO
226    status=NF90_PUT_VAR(ncid,hybmid, val(:))
227 
228    status=NF90_PUT_VAR(ncid,P0id, Preff)
229   
230    status=NF90_CLOSE(ncid)
231  END IF
232  END SUBROUTINE write_apbp
233
234END MODULE disvert_mod
Note: See TracBrowser for help on using the repository browser.