source: codes/icosagcm/trunk/src/vertical/disvert.f90 @ 1057

Last change on this file since 1057 was 1057, checked in by millour, 3 years ago

Add a "plugin" option to disvert. Only the interface to the plugin is present inthe realted disvert_plugin.f90; a placeholder for a user-provided routine set byan external driver.
EM
File size: 9.5 KB
Line 
1MODULE disvert_mod
2  USE icosa
3  USE abort_mod
4  REAL(rstd), SAVE, POINTER :: ap(:)
5!$OMP THREADPRIVATE(ap)
6  REAL(rstd), SAVE, POINTER :: bp(:)
7!$OMP THREADPRIVATE(bp)
8  REAL(rstd), SAVE, POINTER :: presnivs(:)
9!$OMP THREADPRIVATE(presnivs)
10  REAL(rstd), SAVE, POINTER :: mass_al(:), mass_bl(:), mass_ak(:), mass_bk(:), mass_dak(:), mass_dbk(:)
11!$OMP THREADPRIVATE(mass_al, mass_bl, mass_ak, mass_bk, mass_dak, mass_dbk)
12
13  REAL(rstd) :: ptop ! pressure at top of atmosphere l=llm+1
14!$OMP THREADPRIVATE(ptop)
15  ! vertical coordinate : Lagrangian or mass-based
16  INTEGER, SAVE :: caldyn_eta
17!$OMP THREADPRIVATE(caldyn_eta)
18  INTEGER, PARAMETER :: eta_mass=1, eta_lag=2
19  LOGICAL,SAVE :: ap_bp_present
20!$OMP THREADPRIVATE(ap_bp_present)
21
22CONTAINS
23
24  SUBROUTINE init_disvert
25  USE disvert_std_mod, ONLY: ap_std=>ap, bp_std=>bp, presnivs_std=>presnivs, init_disvert_std=>init_disvert
26  USE disvert_strato_mod, ONLY: ap_strato=>ap, bp_strato=>bp, presnivs_strato=>presnivs, init_disvert_strato=>init_disvert, init_disvert_strato_custom
27  USE disvert_apbp_mod, ONLY: ap_apbp=>ap, bp_apbp=>bp, presnivs_apbp=>presnivs, init_disvert_apbp=>init_disvert
28  USE disvert_ncar_mod, ONLY: ap_ncar=>ap, bp_ncar=>bp, presnivs_ncar=>presnivs, init_disvert_ncar=>init_disvert
29  USE disvert_ncarl30_mod, ONLY: ap_ncarl30=>ap, bp_ncarl30=>bp, presnivs_ncarl30=>presnivs, init_disvert_ncarl30=>init_disvert
30  USE disvert_dcmip31_mod, ONLY: ap_dcmip31=>ap, bp_dcmip31=>bp, presnivs_dcmip31=>presnivs, init_disvert_dcmip31=>init_disvert
31  USE disvert_dcmip200_mod, ONLY: ap_dcmip200=>ap, bp_dcmip200=>bp, presnivs_dcmip200=>presnivs, init_disvert_dcmip200=>init_disvert
32  USE disvert_plugin_mod, ONLY: ap_plugin=>ap, bp_plugin=>bp, presnivs_plugin=>presnivs, init_disvert_plugin=>init_disvert
33  USE icosa
34  USE mpipara
35  IMPLICIT NONE
36    CHARACTER(len=255) :: def
37    INTEGER :: l
38
39    def='eta_mass'
40    CALL getin('caldyn_eta',def)
41    SELECT CASE(TRIM(def))
42    CASE('eta_mass')
43       caldyn_eta=eta_mass
44    CASE('eta_lag')
45       caldyn_eta=eta_lag
46    CASE DEFAULT
47       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_eta : <', TRIM(def),'> options are <eta_mass>, <eta_lag>'
48       STOP
49    END SELECT
50    IF (is_mpi_root) PRINT *, 'caldyn_eta=',TRIM(def), caldyn_eta, eta_lag, eta_mass
51
52    def='std'
53    ap_bp_present=.TRUE.
54    CALL getin("disvert",def)
55
56    SELECT CASE (TRIM(def))
57      CASE('none')
58         IF(caldyn_eta/=eta_lag) THEN
59            IF(is_mpi_root) PRINT *,'disvert_type=none must be used with caldyn_eta=eta_lag'
60            STOP
61         END IF
62         ap_bp_present=.FALSE.
63      CASE('std')
64   
65        CALL init_disvert_std
66        ap=>ap_std
67        bp=>bp_std
68        presnivs=>presnivs_std
69     
70      CASE('strato')
71       
72        CALL init_disvert_strato
73        ap=>ap_strato
74        bp=>bp_strato
75        presnivs=>presnivs_strato
76
77      CASE('strato_custom')
78       
79        CALL init_disvert_strato_custom
80        ap=>ap_strato
81        bp=>bp_strato
82        presnivs=>presnivs_strato
83
84      CASE('read_apbp')
85   
86        CALL init_disvert_apbp
87        ap=>ap_apbp
88        bp=>bp_apbp
89        presnivs=>presnivs_apbp
90       
91      CASE('plugin')
92       
93        CALL init_disvert_plugin
94        ap=>ap_plugin
95        bp=>bp_plugin
96        presnivs=>presnivs_plugin
97       
98      CASE ('ncar')
99
100        CALL init_disvert_ncar
101        ap=>ap_ncar
102        bp=>bp_ncar
103        presnivs=>presnivs_ncar
104
105      CASE ('dcmip31')
106
107        CALL init_disvert_dcmip31
108        ap=>ap_dcmip31
109        bp=>bp_dcmip31
110        presnivs=>presnivs_dcmip31
111
112      CASE ('dcmip200')
113
114        CALL init_disvert_dcmip200
115        ap=>ap_dcmip200
116        bp=>bp_dcmip200
117        presnivs=>presnivs_dcmip200
118
119      CASE ('ncarl30')
120
121        CALL init_disvert_ncarl30
122        ap=>ap_ncarl30
123        bp=>bp_ncarl30
124        presnivs=>presnivs_ncarl30
125      CASE default
126        IF (is_mpi_root) PRINT*,'Bad selector for variable disvert : <', TRIM(def),"> options are <std>, <ncar>, <ncarl30>" 
127        STOP
128       
129    END SELECT
130
131    IF(ap_bp_present) THEN
132       ! Convert from pressure coordinate to mass coordinate
133       ! pk = ap(k) + bp(k)*ps(ij) = ptop + g*( mass_ak(k) + mass_bk(k)*ms(ij) )
134       ptop = ap(llm+1)
135       ALLOCATE(mass_al(llm+1))
136       ALLOCATE(mass_bl(llm+1))
137       ALLOCATE(mass_ak(llm))
138       ALLOCATE(mass_bk(llm))
139       ALLOCATE(mass_dak(llm))
140       ALLOCATE(mass_dbk(llm))
141
142       ! FIXME : leave ps for the moment ; change ps to ms later
143       DO l=1,llm+1
144          !       mass_al(l) = (ap(l)-ptop)/g
145          !       mass_bl(l) = bp(l)/g
146          mass_al(l) = (ap(l)-ptop)
147          mass_bl(l) = bp(l)
148       END DO
149       DO l=1,llm
150          mass_ak(l) = .5*(mass_al(l)+mass_al(l+1))
151          mass_bk(l) = .5*(mass_bl(l)+mass_bl(l+1))
152          mass_dak(l) = mass_al(l)-mass_al(l+1)  ! >0
153          mass_dbk(l) = mass_bl(l)-mass_bl(l+1)  ! >0
154       END DO
155    END IF
156
157    !$acc enter data copyin(ap(:), bp(:), mass_dak(:), mass_dbk(:)) async
158   
159  END SUBROUTINE init_disvert 
160 
161  SUBROUTINE compute_rhodz(comp, ps, rhodz, ondevice)
162    USE icosa
163    USE omp_para
164    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
165    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
166    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
167    LOGICAL, INTENT(IN), OPTIONAL :: ondevice ! .TRUE. compute on device, .FALSE. stay on host
168    REAL(rstd) :: m, err
169    INTEGER :: l,i,j,ij,dd
170    LOGICAL :: ondevice_
171
172    err=0.
173
174    ! by default, do not compute on device
175    ondevice_ = .FALSE.
176    IF (PRESENT(ondevice)) ondevice_ = ondevice
177
178    IF(comp) THEN
179       dd=1
180    ELSE
181       dd=0
182    END IF
183
184    !$acc data present(ps(:), rhodz(:,:), ap(:), bp(:)) async if (ondevice_)
185
186    IF(comp) THEN
187      !$acc parallel loop collapse(3) async if (ondevice_)
188      DO l = ll_begin, ll_end
189        DO j=jj_begin-dd,jj_end+dd
190          DO i=ii_begin-dd,ii_end+dd
191            ij=(j-1)*iim+i
192            m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
193            rhodz(ij,l) = m
194          ENDDO
195        ENDDO
196      ENDDO
197    ELSE
198      !$acc parallel loop reduction(max:err) collapse(3) async if (ondevice_)
199      DO l = ll_begin, ll_end
200        DO j=jj_begin-dd,jj_end+dd
201          DO i=ii_begin-dd,ii_end+dd
202            ij=(j-1)*iim+i
203            m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
204            err = MAX(err,abs(m-rhodz(ij,l)))
205          ENDDO
206        ENDDO
207      ENDDO
208
209      !$acc wait
210      IF(err>1e-10) THEN
211          PRINT *, 'Discrepancy between ps and rhodz detected', err
212          STOP
213      ENDIF
214    ENDIF
215
216    !$acc end data
217  END SUBROUTINE compute_rhodz
218
219
220  SUBROUTINE write_apbp
221  USE icosa
222  USE netcdf_mod
223  IMPLICIT NONE
224    REAL(rstd) :: val(llm)
225    INTEGER :: status
226    INTEGER :: lev,ilev
227    INTEGER :: ncid,levid,ilevid,hyaiid,hybiid,hyamid,hybmid,P0id
228    INTEGER :: l
229   
230    IF (no_io) RETURN
231   
232!$OMP BARRIER
233!$OMP MASTER   
234    status = NF90_CREATE('apbp.nc', NF90_CLOBBER, ncid)
235    status = NF90_DEF_DIM(ncid,'lev',llm,lev)
236    status = NF90_DEF_DIM(ncid,'ilev',llm+1,ilev)
237   
238    status = NF90_DEF_VAR(ncid,'lev',NF90_DOUBLE,(/ lev /),levid)
239    status = NF90_PUT_ATT(ncid,levid,"long_name","hybrid level at midpoints (1000*(A+B))")
240    status = NF90_PUT_ATT(ncid,levid,"units","Pa")
241    status = NF90_PUT_ATT(ncid,levid,"positive","up")
242    status = NF90_PUT_ATT(ncid,levid,"standard_name","atmosphere_hybrid_sigma_pressure_coordinate")
243    status = NF90_PUT_ATT(ncid,levid,"formula_terms","a: hyam b: hybm p0: P0 ps: PS")
244
245    status = NF90_DEF_VAR(ncid,'ilev',NF90_DOUBLE,(/ ilev /),ilevid)
246    status = NF90_PUT_ATT(ncid,ilevid,"long_name","hybrid level at interfaces (1000*(A+B))")
247    status = NF90_PUT_ATT(ncid,ilevid,"units","Pa")
248    status = NF90_PUT_ATT(ncid,ilevid,"positive","up")
249    status = NF90_PUT_ATT(ncid,ilevid,"standard_name","atmosphere_hybrid_sigma_pressure_coordinate")
250    status = NF90_PUT_ATT(ncid,ilevid,"formula_terms","a: hyai b: hybi p0: P0 ps: PS")
251   
252    status = NF90_DEF_VAR(ncid,'hyai',NF90_DOUBLE,(/ ilev /),hyaiid)
253    status = NF90_PUT_ATT(ncid,hyaiid,"long_name","hybrid A coefficient at layer interfaces")
254
255    status = NF90_DEF_VAR(ncid,'hybi',NF90_DOUBLE,(/ ilev /),hybiid)
256    status = NF90_PUT_ATT(ncid,hybiid,"long_name","hybrid B coefficient at layer interfaces")
257
258    status = NF90_DEF_VAR(ncid,'hyam',NF90_DOUBLE,(/ lev /),hyamid)
259    status = NF90_PUT_ATT(ncid,hyamid,"long_name","hybrid A coefficient at midpoint interfaces")
260
261    status = NF90_DEF_VAR(ncid,'hybm',NF90_DOUBLE,(/ lev /),hybmid)
262    status = NF90_PUT_ATT(ncid,hybmid,"long_name","hybrid B coefficient at midpoint interfaces")
263   
264    status = NF90_DEF_VAR(ncid,'P0',NF90_DOUBLE,varid=P0id)
265
266    status = NF90_ENDDEF(ncid)   
267
268    IF(ap_bp_present) THEN
269   
270    status=NF90_PUT_VAR(ncid,ilevid, ap(:)+bp(:)*Preff)
271   
272    DO l=1,llm
273      val(l)= 0.5*(ap(l+1)+ap(l)+Preff*(bp(l)+bp(l+1)))
274    ENDDO
275   
276    status=NF90_PUT_VAR(ncid,levid, val)
277
278    status=NF90_PUT_VAR(ncid,hyaiid, ap(:)/Preff)
279    status=NF90_PUT_VAR(ncid,hybiid, bp(:))
280   
281    DO l=1,llm
282      val(l)= 0.5*(ap(l+1)+ap(l))
283    ENDDO
284    status=NF90_PUT_VAR(ncid,hyamid, val(:)/Preff)
285   
286    DO l=1,llm
287      val(l)= 0.5*(bp(l+1)+bp(l))
288    ENDDO
289    status=NF90_PUT_VAR(ncid,hybmid, val(:))
290 
291    status=NF90_PUT_VAR(ncid,P0id, Preff)
292   
293    END IF
294
295    status=NF90_CLOSE(ncid)
296!$OMP END MASTER
297!$OMP BARRIER
298  END SUBROUTINE write_apbp
299
300END MODULE disvert_mod
Note: See TracBrowser for help on using the repository browser.