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

Last change on this file since 953 was 953, checked in by adurocher, 5 years ago

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