source: LMDZ6/branches/contrails/libf/phylmdiso/phyaqua_mod.F90 @ 5440

Last change on this file since 5440 was 5285, checked in by abarral, 8 weeks ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules
File size: 30.5 KB
Line 
1!
2! $Id: phyaqua_mod.F90 3579 2019-10-09 13:11:07Z fairhead $
3!
4MODULE phyaqua_mod
5  ! Routines complementaires pour la physique planetaire.
6  IMPLICIT NONE
7
8CONTAINS
9
10  SUBROUTINE iniaqua(nlon,year_len,iflag_phys)
11
12    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13    ! Creation d'un etat initial et de conditions aux limites
14    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
15    ! du modele LMDZ dans sa version terrestre.
16    ! iflag_phys est un parametre qui controle
17    ! iflag_phys = N
18    ! de 100 a 199 : aqua planetes avec SST forcees
19    ! N-100 determine le type de SSTs
20    ! de 200 a 299 : terra planetes avec Ts calcule
21
22    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23
24    USE dimphy, ONLY: klon
25    USE geometry_mod, ONLY : latitude
26    USE surface_data, ONLY: type_ocean, ok_veget
27    USE pbl_surface_mod, ONLY: pbl_surface_init
28    USE fonte_neige_mod, ONLY: fonte_neige_init
29#ifdef ISO
30  USE fonte_neige_mod,  ONLY : fonte_neige_init_iso
31  USE pbl_surface_mod,  ONLY : pbl_surface_init_iso
32#endif
33    USE phys_state_var_mod
34    USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
35                                day_ini,day_end
36    USE indice_sol_mod
37    USE nrtype, ONLY: pi
38!    USE ioipsl
39    USE mod_phys_lmdz_para, ONLY: is_master
40    USE mod_phys_lmdz_transfert_para, ONLY: bcast
41    USE mod_grid_phy_lmdz
42    USE ioipsl_getin_p_mod, ONLY : getin_p
43    USE phys_cal_mod , ONLY: calend, year_len_phy => year_len
44#ifdef ISO
45  USE infotrac_phy, ONLY: niso
46#endif
47    USE clesphys_mod_h
48    USE dimsoil_mod_h, ONLY: nsoilmx
49
50    USE yomcst_mod_h
51IMPLICIT NONE
52
53
54
55    INTEGER, INTENT (IN) :: nlon, year_len, iflag_phys
56    ! IM ajout latfi, lonfi
57!    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
58
59    INTEGER type_profil, type_aqua
60
61    ! Ajouts initialisation des surfaces
62    REAL :: run_off_lic_0(nlon)
63    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
64    REAL :: tsoil(nlon, nsoilmx, nbsrf)
65    REAL :: tslab(nlon), seaice(nlon)
66    REAL fder(nlon)
67
68#ifdef ISO
69    REAL :: xtrun_off_lic_0(niso,nlon)
70    REAL :: xtsolsrf(niso,nlon, nbsrf), xtsnsrf(niso,nlon, nbsrf)
71    REAL :: Rland_ice(niso,nlon)
72#endif
73
74
75
76    ! Arguments :
77    ! -----------
78
79    ! integer radpas
80    INTEGER it, unit, i, k, itap
81
82    REAL rugos, albedo
83    REAL tsurf
84    REAL time, timestep, day, day0
85    REAL qsol_f
86    REAL rugsrel(nlon)
87    LOGICAL alb_ocean
88
89    CHARACTER *80 ans, file_forctl, file_fordat, file_start
90    CHARACTER *100 file, var
91    CHARACTER *2 cnbl
92
93    REAL phy_nat(nlon, year_len)
94    REAL phy_alb(nlon, year_len)
95    REAL phy_sst(nlon, year_len)
96    REAL phy_bil(nlon, year_len)
97    REAL phy_rug(nlon, year_len)
98    REAL phy_ice(nlon, year_len)
99    REAL phy_fter(nlon, year_len)
100    REAL phy_foce(nlon, year_len)
101    REAL phy_fsic(nlon, year_len)
102    REAL phy_flic(nlon, year_len)
103
104    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
105!$OMP THREADPRIVATE(read_climoz)
106
107    ! -------------------------------------------------------------------------
108    ! declaration pour l'appel a phyredem
109    ! -------------------------------------------------------------------------
110
111    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
112    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
113    ! real pbl_tke(nlon,llm,nbsrf)
114    ! real rain_fall(nlon),snow_fall(nlon)
115    ! real solsw(nlon), sollw(nlon),radsol(nlon)
116    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
117    ! real ratqs(nlon,llm)
118    ! real clwcon(nlon,llm)
119
120    INTEGER longcles
121    PARAMETER (longcles=20)
122    REAL clesphy0(longcles)
123
124
125    ! -----------------------------------------------------------------------
126    ! dynamial tendencies :
127    ! ---------------------
128
129    INTEGER l, ierr, aslun
130
131    REAL paire
132
133    ! Local
134    CHARACTER (LEN=20) :: modname='phyaqua'
135    CHARACTER (LEN=80) :: abort_message
136
137
138    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
139    ! INITIALISATIONS
140    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
141
142    ! -----------------------------------------------------------------------
143    ! Initialisations  des constantes
144    ! -------------------------------
145
146    !IF (calend .EQ. "earth_360d") Then
147      year_len_phy = year_len
148    !END IF
149   
150    if (year_len.ne.360) then
151      write (*,*) year_len
152      write (*,*) 'iniaqua: 360 day calendar is required !'
153      stop
154    endif
155
156    type_aqua = iflag_phys/100
157    type_profil = iflag_phys - type_aqua*100
158    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
159
160    IF (klon/=nlon) THEN
161      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
162      abort_message= 'probleme de dimensions dans iniaqua'
163      CALL abort_physic(modname,abort_message,1)
164    END IF
165    CALL phys_state_var_init(read_climoz)
166
167
168    read_climoz = 0
169    day0 = 217.
170    day = day0
171    it = 0
172    time = 0.
173
174    ! -----------------------------------------------------------------------
175    ! initialisations de la physique
176    ! -----------------------------------------------------------------------
177
178    day_ini = day_ref
179    day_end = day_ini + ndays
180
181    nbapp_rad = 24
182    CALL getin_p('nbapp_rad', nbapp_rad)
183
184    ! ---------------------------------------------------------------------
185    ! Creation des conditions aux limites:
186    ! ------------------------------------
187    ! Initialisations des constantes
188    ! Ajouter les manquants dans planete.def... (albedo etc)
189    co2_ppm = 348.
190    CALL getin_p('co2_ppm', co2_ppm)
191
192    solaire = 1365.
193    CALL getin_p('solaire', solaire)
194 
195    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
196    ! type_aqua
197    alb_ocean = .TRUE.
198    CALL getin_p('alb_ocean', alb_ocean)
199
200    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
201    WRITE (*, *) 'iniaqua: solaire=', solaire
202    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
203
204    radsol = 0.
205    qsol_f = 10.
206
207    ! Conditions aux limites:
208    ! -----------------------
209
210    qsol(:) = qsol_f
211    rugsrel = 0.0 ! (rugsrel = rugoro)
212    rugoro = 0.0
213    u_ancien = 0.0
214    v_ancien = 0.0
215    agesno = 50.0
216    ! Relief plat
217    zmea = 0.
218    zstd = 0.
219    zsig = 0.
220    zgam = 0.
221    zthe = 0.
222    zpic = 0.
223    zval = 0.
224
225    ! Une seule surface
226    pctsrf = 0.
227    IF (type_aqua==1) THEN
228      rugos = 1.E-4
229      albedo = 0.19
230      pctsrf(:, is_oce) = 1.
231    ELSE IF (type_aqua==2) THEN
232      rugos = 0.03
233      albedo = 0.1
234      pctsrf(:, is_ter) = 1.
235    END IF
236
237    CALL getin_p('rugos', rugos)
238
239    WRITE (*, *) 'iniaqua: rugos=', rugos
240    zmasq(:) = pctsrf(:, is_ter)
241
242    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
243    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
244
245    ! Si alb_ocean on calcule un albedo oceanique moyen
246    ! if (alb_ocean) then
247    ! Voir pourquoi on avait ca.
248    ! CALL ini_alb_oce(phy_alb)
249    ! else
250    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
251    ! endif !alb_ocean
252
253    DO i = 1, year_len
254      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
255      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
256
257      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
258      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
259      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
260      phy_ice(:, i) = 0.0 ! fraction de glace (?)
261      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
262      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
263      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
264      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
265    END DO
266    ! IM calcul profil sst
267    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
268
269    IF (grid_type==unstructured) THEN
270      CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
271                             phy_fter, phy_foce, phy_flic, phy_fsic)
272    ELSE
273     
274       CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
275                     phy_fter, phy_foce, phy_flic, phy_fsic)
276    ENDIF
277
278    ! ---------------------------------------------------------------------
279    ! Ecriture de l'etat initial:
280    ! ---------------------------
281
282
283    ! Ecriture etat initial physique
284
285    timestep = pdtphys
286    radpas = nint(rday/timestep/float(nbapp_rad))
287
288    DO i = 1, longcles
289      clesphy0(i) = 0.
290    END DO
291    clesphy0(1) = float(iflag_con)
292    clesphy0(2) = float(nbapp_rad)
293    ! IF( cycle_diurne  ) clesphy0(3) =  1.
294    clesphy0(3) = 1. ! cycle_diurne
295    clesphy0(4) = 1. ! soil_model
296    clesphy0(5) = 1. ! new_oliq
297    clesphy0(6) = 0. ! ok_orodr
298    clesphy0(7) = 0. ! ok_orolf
299    clesphy0(8) = 0. ! ok_limitvrai
300
301
302    ! =======================================================================
303    ! Profils initiaux
304    ! =======================================================================
305
306    ! On initialise les temperatures de surfaces comme les sst
307    DO i = 1, nlon
308      ftsol(i, :) = phy_sst(i, 1)
309      tsoil(i, :, :) = phy_sst(i, 1)
310      tslab(i) = phy_sst(i, 1)
311    END DO
312
313    falbe(:, :) = albedo
314    falblw(:, :) = albedo
315    rain_fall(:) = 0.
316    snow_fall(:) = 0.
317    solsw(:) = 0.
318    sollw(:) = 0.
319    radsol(:) = 0.
320
321    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
322    ! intialisation bidon mais pas grave
323    t_ancien(:, :) = 0.
324    q_ancien(:, :) = 0.
325    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
326    rnebcon = 0.
327    ratqs = 0.
328    clwcon = 0.
329    pbl_tke = 1.E-8
330
331    ! variables supplementaires pour appel a plb_surface_init
332    fder(:) = 0.
333    seaice(:) = 0.
334    run_off_lic_0 = 0.
335    fevap = 0.
336#ifdef ISO
337    xtrun_off_lic_0 = 0.
338#endif
339
340
341    ! Initialisations necessaires avant phyredem
342    type_ocean = 'force'
343    CALL fonte_neige_init(run_off_lic_0)
344#ifdef ISO
345   CALL fonte_neige_init_iso(xtrun_off_lic_0)
346#endif
347    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
348    snsrf(:, :) = 0. ! couverture de neige des sous surface
349    z0m(:, :) = rugos ! couverture de neige des sous surface
350    z0h=z0m
351
352
353    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
354#ifdef ISO
355  CALL pbl_surface_init_iso(xtsnsrf,Rland_ice)
356#endif
357
358    PRINT *, 'iniaqua: before phyredem'
359
360    pbl_tke(:,:,:) = 1.e-8
361    falb1 = albedo
362    falb2 = albedo
363    zmax0 = 0.
364    f0 = 0.
365    sig1 = 0.
366    w01 = 0.
367    wake_deltat = 0.
368    wake_deltaq = 0.
369    wake_s = 0.
370    wake_dens = 0.
371    wake_cstar = 0.
372    wake_pe = 0.
373    wake_fip = 0.
374    fm_therm = 0.
375    entr_therm = 0.
376    detr_therm = 0.
377    ale_bl = 0.
378    ale_bl_trig =0.
379    alp_bl =0.
380    treedrg(:,:,:)=0.
381
382    u10m = 0.
383    v10m = 0.
384
385    ql_ancien   = 0.
386    qs_ancien   = 0.
387    u_ancien    = 0.
388    v_ancien    = 0.
389    prw_ancien  = 0.
390    prlw_ancien = 0.
391    prsw_ancien = 0. 
392
393    ale_wake    = 0.
394    ale_bl_stat = 0. 
395
396
397!ym error : the sub surface dimension is the third not second : forgotten for iniaqua
398!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
399!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
400    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
401    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6
402
403!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
404!ym probably the uninitialized value was 0 for standard (regular grid) case
405    falb_dif(:,:,:)=0
406
407
408    CALL phyredem('startphy.nc')
409
410    PRINT *, 'iniaqua: after phyredem'
411    CALL phys_state_var_end
412
413    RETURN
414  END SUBROUTINE iniaqua
415
416
417  ! ====================================================================
418  ! ====================================================================
419  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
420    USE yomcst_mod_h
421    USE dimphy
422    IMPLICIT NONE
423    ! ====================================================================
424    ! =============================================================
425    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
426    ! Auteur : A. Campoy et F. Hourdin
427    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
428    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
429    ! connaissant la declinaison, la latitude et la longitude.
430
431    ! Dans cette version particuliere, on calcule le rayonnement
432    ! moyen sur l'année à chaque latitude.
433    ! angle zenithal calculé pour obtenir un
434    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
435    ! en moyenne annuelle.
436    ! Spécifique de la terre. Utilisé pour les aqua planetes.
437
438    ! Rque   : Different de la routine angle en ce sens que zenang
439    ! fournit des moyennes de pmu0 et non des valeurs
440    ! instantanees, du coup frac prend toutes les valeurs
441    ! entre 0 et 1.
442    ! Date   : premiere version le 13 decembre 1994
443    ! revu pour  GCM  le 30 septembre 1996
444    ! ===============================================================
445    ! longi----INPUT : la longitude vraie de la terre dans son plan
446    ! solaire a partir de l'equinoxe de printemps (degre)
447    ! gmtime---INPUT : temps universel en fraction de jour
448    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
449    ! lat------INPUT : latitude en degres
450    ! long-----INPUT : longitude en degres
451    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
452    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
453    ! ================================================================
454    ! ================================================================
455    LOGICAL cycle_diurne
456    REAL gmtime
457    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
458    ! ================================================================
459    INTEGER i
460    REAL gmtime1, gmtime2
461    REAL pi_local
462
463
464    REAL rmu0m(klon), rmu0a(klon)
465
466
467    pi_local = 4.0*atan(1.0)
468
469    ! ================================================================
470    ! Calcul de l'angle zenithal moyen sur la journee
471    ! ================================================================
472
473    DO i = 1, klon
474      fract(i) = 1.
475      ! Calcule du flux moyen
476      IF (abs(rlat(i))<=28.75) THEN
477        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
478      ELSE IF (abs(rlat(i))<=43.75) THEN
479        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
480      ELSE IF (abs(rlat(i))<=71.25) THEN
481        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
482      ELSE
483        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
484      END IF
485    END DO
486
487    ! ================================================================
488    ! Avec ou sans cycle diurne
489    ! ================================================================
490
491    IF (cycle_diurne) THEN
492
493      ! On redecompose flux  au sommet suivant un cycle diurne idealise
494      ! identique a toutes les latitudes.
495
496      DO i = 1, klon
497        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
498        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
499          rmu0a(i)/sqrt(2.)
500      END DO
501
502      DO i = 1, klon
503        IF (rmu0(i)<=0.) THEN
504          rmu0(i) = 0.
505          fract(i) = 0.
506        ELSE
507          fract(i) = 1.
508        END IF
509      END DO
510
511      ! Affichage de l'angel zenitale
512      ! print*,'************************************'
513      ! print*,'************************************'
514      ! print*,'************************************'
515      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
516      ! print*,'rmu0m=',rmu0m(i)
517      ! print*,'rmu0a=',rmu0a(i)
518      ! print*,'rmu0=',rmu0(i)
519
520    ELSE
521
522      DO i = 1, klon
523        fract(i) = 0.5
524        rmu0(i) = rmu0m(i)*2.
525      END DO
526
527    END IF
528
529    RETURN
530  END SUBROUTINE zenang_an
531
532  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
533
534  SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
535      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
536
537    USE mod_phys_lmdz_para, ONLY: is_omp_master, klon_mpi
538    USE mod_phys_lmdz_transfert_para, ONLY: gather_omp
539    USE lmdz_xios
540    IMPLICIT NONE
541
542    INTEGER, INTENT (IN) :: klon
543    REAL, INTENT (IN) :: phy_nat(klon, 360)
544    REAL, INTENT (IN) :: phy_alb(klon, 360)
545    REAL, INTENT (IN) :: phy_sst(klon, 360)
546    REAL, INTENT (IN) :: phy_bil(klon, 360)
547    REAL, INTENT (IN) :: phy_rug(klon, 360)
548    REAL, INTENT (IN) :: phy_ice(klon, 360)
549    REAL, INTENT (IN) :: phy_fter(klon, 360)
550    REAL, INTENT (IN) :: phy_foce(klon, 360)
551    REAL, INTENT (IN) :: phy_flic(klon, 360)
552    REAL, INTENT (IN) :: phy_fsic(klon, 360)
553
554    REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
555      ! on the whole physics grid
556 
557    PRINT *, 'writelim: Ecriture du fichier limit'
558
559    CALL gather_omp(phy_foce, phy_mpi)
560    IF (is_omp_master) CALL xios_send_field('foce_limout',phy_mpi)
561
562    CALL gather_omp(phy_fsic, phy_mpi)
563    IF (is_omp_master) CALL xios_send_field('fsic_limout',phy_mpi)
564     
565    CALL gather_omp(phy_fter, phy_mpi)
566    IF (is_omp_master) CALL xios_send_field('fter_limout',phy_mpi)
567     
568    CALL gather_omp(phy_flic, phy_mpi)
569    IF (is_omp_master) CALL xios_send_field('flic_limout',phy_mpi)
570
571    CALL gather_omp(phy_sst, phy_mpi)
572    IF (is_omp_master) CALL xios_send_field('sst_limout',phy_mpi)
573
574    CALL gather_omp(phy_bil, phy_mpi)
575    IF (is_omp_master) CALL xios_send_field('bils_limout',phy_mpi)
576
577    CALL gather_omp(phy_alb, phy_mpi)
578    IF (is_omp_master) CALL xios_send_field('alb_limout',phy_mpi)
579
580    CALL gather_omp(phy_rug, phy_mpi)
581    IF (is_omp_master) CALL xios_send_field('rug_limout',phy_mpi)
582
583  END SUBROUTINE writelim_unstruct
584
585
586
587  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
588      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
589
590    USE mod_phys_lmdz_para, ONLY: is_master
591    USE mod_grid_phy_lmdz, ONLY: klon_glo
592    USE mod_phys_lmdz_transfert_para, ONLY: gather
593    USE phys_cal_mod, ONLY: year_len
594    USE netcdf, ONLY: nf90_clobber, nf90_close, nf90_noerr, nf90_strerror, nf90_put_att, nf90_def_var, &
595            nf90_def_dim, nf90_create, nf90_put_var, nf90_unlimited, nf90_global, nf90_64bit_offset, &
596            nf90_enddef
597    USE lmdz_cppkeys_wrapper, ONLY: nf90_format
598    IMPLICIT NONE
599
600    INTEGER, INTENT (IN) :: klon
601    REAL, INTENT (IN) :: phy_nat(klon, year_len)
602    REAL, INTENT (IN) :: phy_alb(klon, year_len)
603    REAL, INTENT (IN) :: phy_sst(klon, year_len)
604    REAL, INTENT (IN) :: phy_bil(klon, year_len)
605    REAL, INTENT (IN) :: phy_rug(klon, year_len)
606    REAL, INTENT (IN) :: phy_ice(klon, year_len)
607    REAL, INTENT (IN) :: phy_fter(klon, year_len)
608    REAL, INTENT (IN) :: phy_foce(klon, year_len)
609    REAL, INTENT (IN) :: phy_flic(klon, year_len)
610    REAL, INTENT (IN) :: phy_fsic(klon, year_len)
611
612    REAL :: phy_glo(klon_glo, year_len) ! temporary variable, to store phy_***(:)
613      ! on the whole physics grid
614    INTEGER :: k
615    INTEGER ierr
616    INTEGER dimfirst(3)
617    INTEGER dimlast(3)
618
619    INTEGER nid, ndim, ntim
620    INTEGER dims(2), debut(2), epais(2)
621    INTEGER id_tim
622    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
623    INTEGER id_fter, id_foce, id_fsic, id_flic
624
625    IF (is_master) THEN
626
627      PRINT *, 'writelim: Ecriture du fichier limit'
628
629      ierr = nf90_create('limit.nc', IOR(nf90_clobber,nf90_64bit_offset), nid)
630
631      ierr = nf90_put_att(nid, nf90_global, 'title', 'Fichier conditions aux limites')
632      ! !        ierr = nf90_def_dim (nid, "points_physiques", klon, ndim)
633      ierr = nf90_def_dim(nid, 'points_physiques', klon_glo, ndim)
634      ierr = nf90_def_dim(nid, 'time', nf90_unlimited, ntim)
635
636      dims(1) = ndim
637      dims(2) = ntim
638
639      ierr = nf90_def_var(nid, 'TEMPS', nf90_format, ntim, id_tim)
640      ierr = nf90_put_att(nid, id_tim, 'title', 'Jour dans l annee')
641
642      ierr = nf90_def_var(nid, 'NAT', nf90_format, dims, id_nat)
643      ierr = nf90_put_att(nid, id_nat, 'title', 'Nature du sol (0,1,2,3)')
644
645      ierr = nf90_def_var(nid, 'SST', nf90_format, dims, id_sst)
646      ierr = nf90_put_att(nid, id_sst, 'title', 'Temperature superficielle de la mer')
647
648      ierr = nf90_def_var(nid, 'BILS', nf90_format, dims, id_bils)
649      ierr = nf90_put_att(nid, id_bils, 'title', 'Reference flux de chaleur au sol')
650
651      ierr = nf90_def_var(nid, 'ALB', nf90_format, dims, id_alb)
652      ierr = nf90_put_att(nid, id_alb, 'title', 'Albedo a la surface')
653
654      ierr = nf90_def_var(nid, 'RUG', nf90_format, dims, id_rug)
655      ierr = nf90_put_att(nid, id_rug, 'title', 'Rugosite')
656
657      ierr = nf90_def_var(nid, 'FTER', nf90_format, dims, id_fter)
658      ierr = nf90_put_att(nid, id_fter, 'title','Frac. Land')
659      ierr = nf90_def_var(nid, 'FOCE', nf90_format, dims, id_foce)
660      ierr = nf90_put_att(nid, id_foce, 'title','Frac. Ocean')
661      ierr = nf90_def_var(nid, 'FSIC', nf90_format, dims, id_fsic)
662      ierr = nf90_put_att(nid, id_fsic, 'title','Frac. Sea Ice')
663      ierr = nf90_def_var(nid, 'FLIC', nf90_format, dims, id_flic)
664      ierr = nf90_put_att(nid, id_flic, 'title','Frac. Land Ice')
665
666      ierr = nf90_enddef(nid)
667      IF (ierr/=nf90_noerr) THEN
668        WRITE (*, *) 'writelim error: failed to end define mode'
669        WRITE (*, *) nf90_strerror(ierr)
670      END IF
671
672
673      ! write the 'times'
674      DO k = 1, year_len
675        ierr = nf90_put_var(nid, id_tim, k, [k])
676        IF (ierr/=nf90_noerr) THEN
677          WRITE (*, *) 'writelim error with temps(k),k=', k
678          WRITE (*, *) nf90_strerror(ierr)
679        END IF
680      END DO
681
682    END IF ! of if (is_master)
683
684    ! write the fields, after having collected them on master
685
686    CALL gather(phy_nat, phy_glo)
687    IF (is_master) THEN
688      ierr = nf90_put_var(nid, id_nat, phy_glo)
689      IF (ierr/=nf90_noerr) THEN
690        WRITE (*, *) 'writelim error with phy_nat'
691        WRITE (*, *) nf90_strerror(ierr)
692      END IF
693    END IF
694
695    CALL gather(phy_sst, phy_glo)
696    IF (is_master) THEN
697      ierr = nf90_put_var(nid, id_sst, phy_glo)
698      IF (ierr/=nf90_noerr) THEN
699        WRITE (*, *) 'writelim error with phy_sst'
700        WRITE (*, *) nf90_strerror(ierr)
701      END IF
702    END IF
703
704    CALL gather(phy_bil, phy_glo)
705    IF (is_master) THEN
706      ierr = nf90_put_var(nid, id_bils, phy_glo)
707      IF (ierr/=nf90_noerr) THEN
708        WRITE (*, *) 'writelim error with phy_bil'
709        WRITE (*, *) nf90_strerror(ierr)
710      END IF
711    END IF
712
713    CALL gather(phy_alb, phy_glo)
714    IF (is_master) THEN
715      ierr = nf90_put_var(nid, id_alb, phy_glo)
716      IF (ierr/=nf90_noerr) THEN
717        WRITE (*, *) 'writelim error with phy_alb'
718        WRITE (*, *) nf90_strerror(ierr)
719      END IF
720    END IF
721
722    CALL gather(phy_rug, phy_glo)
723    IF (is_master) THEN
724      ierr = nf90_put_var(nid, id_rug, phy_glo)
725      IF (ierr/=nf90_noerr) THEN
726        WRITE (*, *) 'writelim error with phy_rug'
727        WRITE (*, *) nf90_strerror(ierr)
728      END IF
729    END IF
730
731    CALL gather(phy_fter, phy_glo)
732    IF (is_master) THEN
733      ierr = nf90_put_var(nid, id_fter, phy_glo)
734      IF (ierr/=nf90_noerr) THEN
735        WRITE (*, *) 'writelim error with phy_fter'
736        WRITE (*, *) nf90_strerror(ierr)
737      END IF
738    END IF
739
740    CALL gather(phy_foce, phy_glo)
741    IF (is_master) THEN
742      ierr = nf90_put_var(nid, id_foce, phy_glo)
743      IF (ierr/=nf90_noerr) THEN
744        WRITE (*, *) 'writelim error with phy_foce'
745        WRITE (*, *) nf90_strerror(ierr)
746      END IF
747    END IF
748
749    CALL gather(phy_fsic, phy_glo)
750    IF (is_master) THEN
751      ierr = nf90_put_var(nid, id_fsic, phy_glo)
752      IF (ierr/=nf90_noerr) THEN
753        WRITE (*, *) 'writelim error with phy_fsic'
754        WRITE (*, *) nf90_strerror(ierr)
755      END IF
756    END IF
757
758    CALL gather(phy_flic, phy_glo)
759    IF (is_master) THEN
760      ierr = nf90_put_var(nid, id_flic, phy_glo)
761      IF (ierr/=nf90_noerr) THEN
762        WRITE (*, *) 'writelim error with phy_flic'
763        WRITE (*, *) nf90_strerror(ierr)
764      END IF
765    END IF
766
767    ! close file:
768    IF (is_master) THEN
769      ierr = nf90_close(nid)
770    END IF
771
772  END SUBROUTINE writelim
773
774  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
775
776  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
777    USE dimphy
778    USE phys_cal_mod , ONLY: year_len
779    IMPLICIT NONE
780
781    INTEGER nlon, type_profil, i, k, j
782    REAL :: rlatd(nlon), phy_sst(nlon, year_len)
783    INTEGER imn, imx, amn, amx, kmn, kmx
784    INTEGER p, pplus, nlat_max
785    PARAMETER (nlat_max=72)
786    REAL x_anom_sst(nlat_max)
787    CHARACTER (LEN=20) :: modname='profil_sst'
788    CHARACTER (LEN=80) :: abort_message
789
790    IF (klon/=nlon) THEN
791       abort_message='probleme de dimensions dans profil_sst'
792       CALL abort_physic(modname,abort_message,1)
793    ENDIF
794    WRITE (*, *) ' profil_sst: type_profil=', type_profil
795    DO i = 1, year_len
796      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
797
798      ! Rajout fbrlmd
799
800      IF (type_profil==1) THEN
801        ! Méthode 1 "Control" faible plateau à l'Equateur
802        DO j = 1, klon
803          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
804          ! PI/3=1.047197551
805          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
806            phy_sst(j, i) = 273.
807          END IF
808        END DO
809      END IF
810      IF (type_profil==2) THEN
811        ! Méthode 2 "Flat" fort plateau à l'Equateur
812        DO j = 1, klon
813          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
814          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
815            phy_sst(j, i) = 273.
816          END IF
817        END DO
818      END IF
819
820
821      IF (type_profil==3) THEN
822        ! Méthode 3 "Qobs" plateau réel à l'Equateur
823        DO j = 1, klon
824          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
825            rlatd(j))**4)
826          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
827            phy_sst(j, i) = 273.
828          END IF
829        END DO
830      END IF
831
832      IF (type_profil==4) THEN
833        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
834        DO j = 1, klon
835          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
836            rlatd(j))**4)
837          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
838            phy_sst(j, i) = 273.
839          END IF
840        END DO
841      END IF
842
843      IF (type_profil==5) THEN
844        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
845        DO j = 1, klon
846          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
847            *rlatd(j))**4)
848          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
849            phy_sst(j, i) = 273. + 2.
850          END IF
851
852        END DO
853      END IF
854
855      IF (type_profil==6) THEN
856        ! Méthode 6 "cst" valeur constante de SST
857        DO j = 1, klon
858          phy_sst(j, i) = 288.
859        END DO
860      END IF
861
862
863      IF (type_profil==7) THEN
864        ! Méthode 7 "cst" valeur constante de SST +2
865        DO j = 1, klon
866          phy_sst(j, i) = 288. + 2.
867        END DO
868      END IF
869
870      p = 0
871      IF (type_profil==8) THEN
872        ! Méthode 8 profil anomalies SST du modèle couplé AR4
873        DO j = 1, klon
874          IF (rlatd(j)==rlatd(j-1)) THEN
875            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
876              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
877            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
878              phy_sst(j, i) = 273. + x_anom_sst(pplus)
879            END IF
880          ELSE
881            p = p + 1
882            pplus = 73 - p
883            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
884              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
885            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
886              phy_sst(j, i) = 273. + x_anom_sst(pplus)
887            END IF
888            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
889          END IF
890        END DO
891      END IF
892
893      IF (type_profil==9) THEN
894        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
895        DO j = 1, klon
896          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
897            *rlatd(j))**4)
898          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
899            phy_sst(j, i) = 273. - 2.
900          END IF
901        END DO
902      END IF
903
904
905      IF (type_profil==10) THEN
906        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
907        DO j = 1, klon
908          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
909            *rlatd(j))**4)
910          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
911            phy_sst(j, i) = 273. + 4.
912          END IF
913        END DO
914      END IF
915
916      IF (type_profil==11) THEN
917        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
918        DO j = 1, klon
919          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
920            rlatd(j))**4)
921          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
922            phy_sst(j, i) = 273.
923          END IF
924        END DO
925      END IF
926
927      IF (type_profil==12) THEN
928        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
929        DO j = 1, klon
930          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
931            *rlatd(j))**4)
932          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
933            phy_sst(j, i) = 273. + 4.
934          END IF
935        END DO
936      END IF
937
938      IF (type_profil==13) THEN
939        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
940        DO j = 1, klon
941          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
942            rlatd(j))**4)
943          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
944            phy_sst(j, i) = 273.
945          END IF
946        END DO
947      END IF
948
949      IF (type_profil==14) THEN
950        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
951        DO j = 1, klon
952          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
953            *rlatd(j))**4)
954          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
955            phy_sst(j, i) = 273.
956          END IF
957        END DO
958      END IF
959
960      if (type_profil.EQ.20) then
961      print*,'Profile SST 20'
962!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
963
964      do j=1,klon
965        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
966      enddo
967      endif
968
969      if (type_profil.EQ.21) then
970      print*,'Profile SST 21'
971!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
972      do j=1,klon
973        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
974      enddo
975      endif
976
977
978
979    END DO
980
981    ! IM beg : verif profil SST: phy_sst
982    amn = min(phy_sst(1,1), 1000.)
983    amx = max(phy_sst(1,1), -1000.)
984    imn = 1
985    kmn = 1
986    imx = 1
987    kmx = 1
988    DO k = 1, year_len
989      DO i = 2, nlon
990        IF (phy_sst(i,k)<amn) THEN
991          amn = phy_sst(i, k)
992          imn = i
993          kmn = k
994        END IF
995        IF (phy_sst(i,k)>amx) THEN
996          amx = phy_sst(i, k)
997          imx = i
998          kmx = k
999        END IF
1000      END DO
1001    END DO
1002
1003    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
1004    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
1005    ! IM end : verif profil SST: phy_sst
1006
1007    RETURN
1008  END SUBROUTINE profil_sst
1009
1010END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.