source: LMDZ6/trunk/libf/phylmd/create_limit_unstruct.F90 @ 4570

Last change on this file since 4570 was 4361, checked in by lguez, 2 years ago

Change calendar attribute "360d" to "360_day"

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