source: LMDZ6/branches/Amaury_dev/libf/phylmd/create_limit_unstruct_mod.F90 @ 5133

Last change on this file since 5133 was 5118, checked in by abarral, 5 months ago

Replace iniprint.h by lmdz_iniprint.f90
(lint) along the way
File size: 11.6 KB
Line 
1MODULE create_limit_unstruct_mod
2  PRIVATE
3  INTEGER, PARAMETER :: lmdep = 12
4
5  PUBLIC create_limit_unstruct
6
7CONTAINS
8
9
10  SUBROUTINE create_limit_unstruct
11    USE dimphy
12    USE lmdz_xios
13    USE ioipsl, ONLY: ioget_year_len
14    USE time_phylmdz_mod, ONLY: annee_ref
15    USE indice_sol_mod
16    USE phys_state_var_mod
17    USE lmdz_phys_para
18    USE lmdz_abort_physic, ONLY: abort_physic
19    USE lmdz_iniprint, ONLY: lunout, prt_level
20    IMPLICIT NONE
21    REAL, DIMENSION(:, :), ALLOCATABLE :: sic
22    REAL, DIMENSION(:, :), ALLOCATABLE :: sst
23    REAL, DIMENSION(klon, lmdep) :: rugos
24    REAL, DIMENSION(klon, lmdep) :: albedo
25    REAL, DIMENSION(:, :), ALLOCATABLE :: sic_mpi
26    REAL, DIMENSION(:, :), ALLOCATABLE :: sst_mpi
27    REAL, DIMENSION(klon_mpi, lmdep) :: rugos_mpi
28    REAL, DIMENSION(klon_mpi, lmdep) :: albedo_mpi
29    INTEGER :: ndays
30    REAL :: fi_ice(klon)
31    REAL, ALLOCATABLE :: sic_year(:, :)
32    REAL, ALLOCATABLE :: sst_year(:, :)
33    REAL, ALLOCATABLE :: rugos_year(:, :)
34    REAL, ALLOCATABLE :: albedo_year(:, :)
35    REAL, ALLOCATABLE :: pctsrf_t(:, :, :)
36    REAL, ALLOCATABLE :: phy_bil(:, :)
37    REAL, ALLOCATABLE :: sst_year_mpi(:, :)
38    REAL, ALLOCATABLE :: rugos_year_mpi(:, :)
39    REAL, ALLOCATABLE :: albedo_year_mpi(:, :)
40    REAL, ALLOCATABLE :: pctsrf_t_mpi(:, :, :)
41    REAL, ALLOCATABLE :: phy_bil_mpi(:, :)
42    INTEGER :: l, k
43    INTEGER :: nbad
44    INTEGER :: sic_time_axis_size
45    INTEGER :: sst_time_axis_size
46    CHARACTER(LEN = 99) :: mess            ! error message
47
48    ndays = ioget_year_len(annee_ref)
49
50    IF (is_omp_master) CALL xios_get_axis_attr("time_sic", n_glo = sic_time_axis_size)
51    CALL bcast_omp(sic_time_axis_size)
52    ALLOCATE(sic_mpi(klon_mpi, sic_time_axis_size))
53    ALLOCATE(sic(klon, sic_time_axis_size))
54
55    IF (is_omp_master) CALL xios_get_axis_attr("time_sst", n_glo = sst_time_axis_size)
56    CALL bcast_omp(sst_time_axis_size)
57    ALLOCATE(sst_mpi(klon_mpi, sst_time_axis_size))
58    ALLOCATE(sst(klon, sst_time_axis_size))
59
60    IF (is_omp_master) THEN
61      CALL xios_recv_field("sic_limit", sic_mpi)
62      CALL xios_recv_field("sst_limit", sst_mpi)
63      CALL xios_recv_field("rugos_limit", rugos_mpi)
64      CALL xios_recv_field("albedo_limit", albedo_mpi)
65    ENDIF
66    CALL scatter_omp(sic_mpi, sic)
67    CALL scatter_omp(sst_mpi, sst)
68    CALL scatter_omp(rugos_mpi, rugos)
69    CALL scatter_omp(albedo_mpi, albedo)
70
71    ALLOCATE(sic_year(klon, ndays))
72    ALLOCATE(sst_year(klon, ndays))
73    ALLOCATE(rugos_year(klon, ndays))
74    ALLOCATE(albedo_year(klon, ndays))
75    ALLOCATE(pctsrf_t(klon, nbsrf, ndays))
76    ALLOCATE(phy_bil(klon, ndays)); phy_bil = 0.0
77
78
79    ! sic
80    IF (sic_time_axis_size==lmdep) THEN
81      CALL time_interpolation(ndays, sic, 'gregorian', sic_year)
82    ELSE IF (sic_time_axis_size==ndays) THEN
83      sic_year = sic
84    ELSE
85      WRITE(mess, *) 'sic time axis is nor montly, nor daily. sic time interpolation ', &
86              'is requiered but is not currently managed'
87      CALL abort_physic('create_limit_unstruct', TRIM(mess), 1)
88    ENDIF
89
90    sic_year(:, :) = sic_year(:, :) / 100.  ! convert percent to fraction
91    WHERE(sic_year(:, :)>1.0) sic_year(:, :) = 1.0    ! Some fractions have some time large negative values
92    WHERE(sic_year(:, :)<0.0) sic_year(:, :) = 0.0    ! probably better to apply alse this filter before horizontal interpolation
93
94    ! sst
95    IF (sst_time_axis_size==lmdep) THEN
96      CALL time_interpolation(ndays, sst, 'gregorian', sst_year)
97    ELSE IF (sst_time_axis_size==ndays) THEN
98      sst_year = sst
99    ELSE
100      WRITE(mess, *)'sic time axis is nor montly, nor daily. sic time interpolation ', &
101              'is requiered but is not currently managed'
102      CALL abort_physic('create_limit_unstruct', TRIM(mess), 1)
103    ENDIF
104    WHERE(sst_year(:, :)<271.38) sst_year(:, :) = 271.38
105
106
107    ! rugos
108    DO l = 1, lmdep
109      WHERE(NINT(zmasq(:))/=1) rugos(:, l) = 0.001
110    ENDDO
111    CALL time_interpolation(ndays, rugos, '360_day', rugos_year)
112
113    ! albedo
114    CALL time_interpolation(ndays, albedo, '360_day', albedo_year)
115
116    DO k = 1, ndays
117      fi_ice = sic_year(:, k)
118      WHERE(fi_ice>=1.0) fi_ice = 1.0
119      WHERE(fi_ice<EPSFRA) fi_ice = 0.0
120      pctsrf_t(:, is_ter, k) = pctsrf(:, is_ter)       ! land soil
121      pctsrf_t(:, is_lic, k) = pctsrf(:, is_lic)       ! land ice
122
123      !!     IF (icefile==trim(fcpldsic)) THEN           ! SIC=pICE*(1-LIC-TER)
124      !!        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
125      !!     ELSE IF (icefile==trim(fhistsic)) THEN      ! SIC=pICE
126      !!        pctsrf_t(:,is_sic,k)=fi_ice(:)
127      !!     ELSE ! icefile==famipsic                    ! SIC=pICE-LIC
128      pctsrf_t(:, is_sic, k) = fi_ice - pctsrf_t(:, is_lic, k)
129      !     END IF
130      WHERE(pctsrf_t(:, is_sic, k)<=0) pctsrf_t(:, is_sic, k) = 0.
131      WHERE(1.0 - zmasq<EPSFRA)
132        pctsrf_t(:, is_sic, k) = 0.0
133        pctsrf_t(:, is_oce, k) = 0.0
134      ELSEWHERE
135        WHERE(pctsrf_t(:, is_sic, k)>=1.0 - zmasq)
136          pctsrf_t(:, is_sic, k) = 1.0 - zmasq
137          pctsrf_t(:, is_oce, k) = 0.0
138        ELSEWHERE
139          pctsrf_t(:, is_oce, k) = 1.0 - zmasq - pctsrf_t(:, is_sic, k)
140          WHERE(pctsrf_t(:, is_oce, k)<EPSFRA)
141            pctsrf_t(:, is_oce, k) = 0.0
142            pctsrf_t(:, is_sic, k) = 1.0 - zmasq
143          END WHERE
144        END WHERE
145      END WHERE
146      nbad = COUNT(pctsrf_t(:, is_oce, k)<0.0)
147      IF(nbad>0) WRITE(lunout, *) 'pb sous maille pour nb point = ', nbad
148      nbad = COUNT(abs(sum(pctsrf_t(:, :, k), dim = 2) - 1.0)>EPSFRA)
149      IF(nbad>0) WRITE(lunout, *) 'pb sous surface pour nb points = ', nbad
150    END DO
151
152    ALLOCATE(sst_year_mpi(klon_mpi, ndays))
153    ALLOCATE(rugos_year_mpi(klon_mpi, ndays))
154    ALLOCATE(albedo_year_mpi(klon_mpi, ndays))
155    ALLOCATE(pctsrf_t_mpi(klon_mpi, nbsrf, ndays))
156    ALLOCATE(phy_bil_mpi(klon_mpi, ndays))
157
158    CALL gather_omp(pctsrf_t, pctsrf_t_mpi)
159    CALL gather_omp(sst_year, sst_year_mpi)
160    CALL gather_omp(phy_bil, phy_bil_mpi)
161    CALL gather_omp(albedo_year, albedo_year_mpi)
162    CALL gather_omp(rugos_year, rugos_year_mpi)
163
164    IF (is_omp_master) THEN
165      CALL xios_send_field("foce_limout", pctsrf_t_mpi(:, is_oce, :))
166      CALL xios_send_field("fsic_limout", pctsrf_t_mpi(:, is_sic, :))
167      CALL xios_send_field("fter_limout", pctsrf_t_mpi(:, is_ter, :))
168      CALL xios_send_field("flic_limout", pctsrf_t_mpi(:, is_lic, :))
169      CALL xios_send_field("sst_limout", sst_year_mpi)
170      CALL xios_send_field("bils_limout", phy_bil_mpi)
171      CALL xios_send_field("alb_limout", albedo_year_mpi)
172      CALL xios_send_field("rug_limout", rugos_year_mpi)
173    ENDIF
174  END SUBROUTINE create_limit_unstruct
175
176
177  SUBROUTINE time_interpolation(ndays, field_in, calendar, field_out)
178    USE lmdz_libmath_pch, ONLY: pchsp_95, pchfe_95
179    USE lmdz_arth, ONLY: arth
180    USE dimphy, ONLY: klon
181    USE ioipsl, ONLY: ioget_year_len
182    USE time_phylmdz_mod, ONLY: annee_ref
183    USE lmdz_phys_para
184    USE lmdz_abort_physic, ONLY: abort_physic
185    USE lmdz_iniprint, ONLY: lunout, prt_level
186    IMPLICIT NONE
187
188    INTEGER, INTENT(IN) :: ndays
189    REAL, INTENT(IN) :: field_in(klon, lmdep)
190    CHARACTER(LEN = *), INTENT(IN) :: calendar
191    REAL, INTENT(OUT) :: field_out(klon, ndays)
192
193    INTEGER :: ndays_in
194    REAL :: timeyear(lmdep)
195    REAL :: yder(lmdep)
196    INTEGER :: ij, ierr, n_extrap
197    LOGICAL :: skip
198
199    CHARACTER (len = 50) :: modname = 'create_limit_unstruct.time_interpolation'
200    CHARACTER (len = 80) :: abort_message
201
202    IF (is_omp_master) ndays_in = year_len(annee_ref, calendar)
203    CALL bcast_omp(ndays_in)
204    IF (is_omp_master) timeyear = mid_months(annee_ref, calendar, lmdep)
205    CALL bcast_omp(timeyear)
206
207    n_extrap = 0
208    skip = .FALSE.
209    DO ij = 1, klon
210      yder = pchsp_95(timeyear, field_in(ij, :), ibeg = 2, iend = 2, vc_beg = 0., vc_end = 0.)
211      CALL pchfe_95(timeyear, field_in(ij, :), yder, skip, arth(0., real(ndays_in) / ndays, ndays), field_out(ij, :), ierr)
212      IF (ierr < 0) THEN
213        abort_message = 'error in pchfe_95'
214        CALL abort_physic(modname, abort_message, 1)
215      endif
216      n_extrap = n_extrap + ierr
217    END DO
218
219    IF (n_extrap /= 0) THEN
220      WRITE(lunout, *) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
221    ENDIF
222
223  END SUBROUTINE time_interpolation
224  !-------------------------------------------------------------------------------
225
226  FUNCTION year_len(y, cal_in)
227
228    !-------------------------------------------------------------------------------
229    USE ioipsl, ONLY: ioget_calendar, ioconf_calendar, lock_calendar, ioget_year_len
230    IMPLICIT NONE
231    !-------------------------------------------------------------------------------
232    ! Arguments:
233    INTEGER :: year_len
234    INTEGER, INTENT(IN) :: y
235    CHARACTER(LEN = *), INTENT(IN) :: cal_in
236    !-------------------------------------------------------------------------------
237    ! Local variables:
238    CHARACTER(LEN = 20) :: cal_out              ! calendar (for outputs)
239    !-------------------------------------------------------------------------------
240    !--- Getting the input calendar to reset at the end of the function
241    CALL ioget_calendar(cal_out)
242
243    !--- Unlocking calendar and setting it to wanted one
244    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
245
246    !--- Getting the number of days in this year
247    year_len = ioget_year_len(y)
248
249    !--- Back to original calendar
250    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
251
252  END FUNCTION year_len
253
254  !-------------------------------------------------------------------------------
255
256
257  !-------------------------------------------------------------------------------
258
259  FUNCTION mid_months(y, cal_in, nm)
260
261    !-------------------------------------------------------------------------------
262    USE ioipsl, ONLY: ioget_calendar, ioconf_calendar, lock_calendar, ioget_mon_len
263    USE lmdz_abort_physic, ONLY: abort_physic
264    IMPLICIT NONE
265    !-------------------------------------------------------------------------------
266    ! Arguments:
267    INTEGER, INTENT(IN) :: y               ! year
268    CHARACTER(LEN = *), INTENT(IN) :: cal_in          ! calendar
269    INTEGER, INTENT(IN) :: nm              ! months/year number
270    REAL, DIMENSION(nm) :: mid_months      ! mid-month times
271    !-------------------------------------------------------------------------------
272    ! Local variables:
273    CHARACTER(LEN = 99) :: mess            ! error message
274    CHARACTER(LEN = 20) :: cal_out         ! calendar (for outputs)
275    INTEGER, DIMENSION(nm) :: mnth            ! months lengths (days)
276    INTEGER :: m               ! months counter
277    INTEGER :: nd              ! number of days
278    INTEGER :: k
279    !-------------------------------------------------------------------------------
280    nd = year_len(y, cal_in)
281
282    IF(nm==12) THEN
283
284      !--- Getting the input calendar to reset at the end of the function
285      CALL ioget_calendar(cal_out)
286
287      !--- Unlocking calendar and setting it to wanted one
288      CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
289
290      !--- Getting the length of each month
291      DO m = 1, nm; mnth(m) = ioget_mon_len(y, m);
292      END DO
293
294      !--- Back to original calendar
295      CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
296
297    ELSE IF(MODULO(nd, nm)/=0) THEN
298      WRITE(mess, '(a,i3,a,i3,a)')'Unconsistent calendar: ', nd, ' days/year, but ', &
299              nm, ' months/year. Months number should divide days number.'
300      CALL abort_physic('mid_months', TRIM(mess), 1)
301
302    ELSE
303      mnth = (/(m, m = 1, nm, nd / nm)/)
304    END IF
305
306    !--- Mid-months times
307    mid_months(1) = 0.5 * REAL(mnth(1))
308    DO k = 2, nm
309      mid_months(k) = mid_months(k - 1) + 0.5 * REAL(mnth(k - 1) + mnth(k))
310    END DO
311
312  END FUNCTION mid_months
313
314
315END MODULE create_limit_unstruct_mod
Note: See TracBrowser for help on using the repository browser.