source: CONFIG/UNIFORM/v7/ICOLMDZOR_v7/SOURCES/DYNAMICO/src/vertical/disvert.f90 @ 5878

Last change on this file since 5878 was 5878, checked in by aclsce, 3 years ago

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